;;; bignum-logical-and bignum-logical-ior bignum-logical-xor
;;; bignum-logical-not bignum-load-byte bignum-deposit-byte
;;; bignum-truncate bignum-plus-p bignum-compare make-small-bignum
-;;; bignum-logcount
+;;; bignum-logbitp bignum-logcount
;;; These symbols define the interface to the compiler:
;;; bignum-type bignum-element-type bignum-index %allocate-bignum
;;; %bignum-length %bignum-set-length %bignum-ref %bignum-set
;;;
;;; PROBLEM 1:
;;; There might be a problem with various LET's and parameters that take a
-;;; digit value. We need to write these so those things stay in 32-bit
+;;; digit value. We need to write these so those things stay in machine
;;; registers and number stack slots. I bind locals to these values, and I
;;; use function on them -- ZEROP, ASH, etc.
;;;
\f
;;;; What's a bignum?
-(eval-when (:compile-toplevel :load-toplevel :execute) ; necessary for DEFTYPE
+(defconstant digit-size sb!vm:n-word-bits)
-(defconstant digit-size 32)
+(defconstant maximum-bignum-length (1- (ash 1 (- sb!vm:n-word-bits
+ sb!vm:n-widetag-bits))))
-(defconstant maximum-bignum-length (1- (ash 1 (- 32 sb!vm:type-bits))))
-
-) ; EVAL-WHEN
+(defconstant all-ones-digit (1- (ash 1 sb!vm:n-word-bits)))
\f
;;;; internal inline routines
(%bignum-length bignum))
;;; %BIGNUM-REF needs to access bignums as obviously as possible, and it needs
-;;; to be able to return 32 bits somewhere no one looks for real objects.
+;;; to be able to return the digit somewhere no one looks for real objects.
(defun %bignum-ref (bignum i)
(declare (type bignum-type bignum)
(type bignum-index i))
(type bignum-index len))
(%digit-0-or-plusp (%bignum-ref bignum (1- len))))
-;;; This should be in assembler, and should not cons intermediate results. It
-;;; returns a 32bit digit and a carry resulting from adding together a, b, and
-;;; an incoming carry.
+;;; This should be in assembler, and should not cons intermediate
+;;; results. It returns a bignum digit and a carry resulting from adding
+;;; together a, b, and an incoming carry.
(defun %add-with-carry (a b carry)
(declare (type bignum-element-type a b)
(type (mod 2) carry))
(%add-with-carry a b carry))
-;;; This should be in assembler, and should not cons intermediate results. It
-;;; returns a 32bit digit and a borrow resulting from subtracting b from a, and
-;;; subtracting a possible incoming borrow.
+;;; This should be in assembler, and should not cons intermediate
+;;; results. It returns a bignum digit and a borrow resulting from
+;;; subtracting b from a, and subtracting a possible incoming borrow.
;;;
;;; We really do: a - b - 1 + borrow, where borrow is either 0 or 1.
(defun %subtract-with-borrow (a b borrow)
(type (mod 2) borrow))
(%subtract-with-borrow a b borrow))
-;;; Multiply two digit-size (32-bit) numbers, returning a 64-bit result
-;;; split into two 32-bit quantities.
+;;; Multiply two digit-size numbers, returning a 2*digit-size result
+;;; split into two digit-size quantities.
(defun %multiply (x y)
(declare (type bignum-element-type x y))
(%multiply x y))
;;; This multiplies x-digit and y-digit, producing high and low digits
;;; manifesting the result. Then it adds the low digit, res-digit, and
-;;; carry-in-digit. Any carries (note, you still have to add two digits at a
-;;; time possibly producing two carries) from adding these three digits get
-;;; added to the high digit from the multiply, producing the next carry digit.
-;;; Res-digit is optional since two uses of this primitive multiplies a single
-;;; digit bignum by a multiple digit bignum, and in this situation there is no
-;;; need for a result buffer accumulating partial results which is where the
-;;; res-digit comes from.
+;;; carry-in-digit. Any carries (note, you still have to add two digits
+;;; at a time possibly producing two carries) from adding these three
+;;; digits get added to the high digit from the multiply, producing the
+;;; next carry digit. Res-digit is optional since two uses of this
+;;; primitive multiplies a single digit bignum by a multiple digit
+;;; bignum, and in this situation there is no need for a result buffer
+;;; accumulating partial results which is where the res-digit comes
+;;; from.
(defun %multiply-and-add (x-digit y-digit carry-in-digit
&optional (res-digit 0))
(declare (type bignum-element-type x-digit y-digit res-digit carry-in-digit))
(declare (type bignum-element-type digit))
(%lognot digit))
-;;; Each of these does the 32-bit unsigned op.
+;;; Each of these does the digit-size unsigned op.
#!-sb-fluid (declaim (inline %logand %logior %logxor))
(defun %logand (a b)
(declare (type bignum-element-type a b))
(declare (type bignum-element-type a b))
(logxor a b))
-;;; This takes a fixnum and sets it up as an unsigned 32-bit quantity. In
-;;; the new system this will mean shifting it right two bits.
+;;; This takes a fixnum and sets it up as an unsigned digit-size
+;;; quantity.
(defun %fixnum-to-digit (x)
(declare (fixnum x))
(logand x (1- (ash 1 digit-size))))
#!-32x16-divide
-;;; This takes three digits and returns the FLOOR'ed result of dividing the
-;;; first two as a 64-bit integer by the third.
+;;; This takes three digits and returns the FLOOR'ed result of
+;;; dividing the first two as a 2*digit-size integer by the third.
;;;
-;;; DO WEIRD let AND setq STUFF TO SLIME THE COMPILER INTO ALLOWING THE %FLOOR
-;;; TRANSFORM TO EXPAND INTO PSEUDO-ASSEMBLER FOR WHICH THE COMPILER CAN LATER
-;;; CORRECTLY ALLOCATE REGISTERS.
+;;; Do weird LET and SETQ stuff to bamboozle the compiler into allowing
+;;; the %FLOOR transform to expand into pseudo-assembler for which the
+;;; compiler can later correctly allocate registers.
(defun %floor (a b c)
(let ((a a) (b b) (c c))
(declare (type bignum-element-type a b c))
;;; unsigned.
(defun %ashr (data count)
(declare (type bignum-element-type data)
- (type (mod 32) count))
+ (type (mod #.sb!vm:n-word-bits) count))
(%ashr data count))
-;;; This takes a 32-bit quantity and shifts it to the left, returning a 32-bit
-;;; quantity.
+;;; This takes a digit-size quantity and shifts it to the left,
+;;; returning a digit-size quantity.
(defun %ashl (data count)
(declare (type bignum-element-type data)
- (type (mod 32) count))
+ (type (mod #.sb!vm:n-word-bits) count))
(%ashl data count))
;;; Do an unsigned (logical) right shift of a digit by Count.
(defun %digit-logical-shift-right (data count)
(declare (type bignum-element-type data)
- (type (mod 32) count))
+ (type (mod #.sb!vm:n-word-bits) count))
(%digit-logical-shift-right data count))
;;; Change the length of bignum to be newlen. Newlen must be the same or
(type bignum-index len))
(%ashr (%bignum-ref bignum (1- len)) (1- digit-size)))
-;;; These take two 32 bit quantities and compare or contrast them without
-;;; wasting time with incorrect type checking.
+;;; These take two digit-size quantities and compare or contrast them
+;;; without wasting time with incorrect type checking.
#!-sb-fluid (declaim (inline %digit-compare %digit-greater))
(defun %digit-compare (x y)
(= x y))
;;; Operations requiring a subtraction without the overhead of intermediate
;;; results, such as GCD, use this. It assumes Result is big enough for the
;;; result.
-(defun subtract-bignum-buffers (a len-a b len-b result)
+(defun subtract-bignum-buffers-with-len (a len-a b len-b result len-res)
(declare (type bignum-type a b)
(type bignum-index len-a len-b))
- (let ((len-res (max len-a len-b)))
- (subtract-bignum-loop a len-a b len-b result len-res
- %normalize-bignum-buffer)))
+ (subtract-bignum-loop a len-a b len-b result len-res
+ %normalize-bignum-buffer))
+
+(defun subtract-bignum-buffers (a len-a b len-b result)
+ (declare (type bignum-type a b)
+ (type bignum-index len-a len-b))
+ (subtract-bignum-loop a len-a b len-b result (max len-a len-b)
+ %normalize-bignum-buffer))
\f
;;;; multiplication
(declare (type bignum-element-type high low))
(if (and (zerop high)
(%digit-0-or-plusp low))
- (let ((low (sb!ext:truly-the (unsigned-byte 31)
+ (let ((low (sb!ext:truly-the (unsigned-byte #.(1- sb!vm:n-word-bits))
(%fixnum-digit-with-correct-sign low))))
(if (eq a-minusp b-minusp)
low
\f
;;;; GCD
-(defun bignum-gcd (a b)
+(eval-when (:compile-toplevel :load-toplevel :execute)
+ ;; The asserts in the GCD implementation are way too expensive to
+ ;; check in normal use, and are disabled here.
+ (sb!xc:defmacro gcd-assert (&rest args)
+ (if nil
+ `(assert ,@args)))
+ ;; We'll be doing a lot of modular arithmetic.
+ (sb!xc:defmacro M (form)
+ `(logand all-ones-digit ,form)))
+
+;;; I'm not sure why I need this FTYPE declaration. Compiled by the
+;;; target compiler, it can deduce the return type fine, but without
+;;; it, we pay a heavy price in BIGNUM-GCD when compiled by the
+;;; cross-compiler. -- CSR, 2004-07-19
+(declaim (ftype (sfunction (bignum-type bignum-index bignum-type bignum-index)
+ sb!vm::positive-fixnum)
+ bignum-factors-of-two))
+(defun bignum-factors-of-two (a len-a b len-b)
+ (declare (type bignum-index len-a len-b) (type bignum-type a b))
+ (do ((i 0 (1+ i))
+ (end (min len-a len-b)))
+ ((= i end) (error "Unexpected zero bignums?"))
+ (declare (type bignum-index i end))
+ (let ((or-digits (%logior (%bignum-ref a i) (%bignum-ref b i))))
+ (unless (zerop or-digits)
+ (return (do ((j 0 (1+ j))
+ (or-digits or-digits (%ashr or-digits 1)))
+ ((oddp or-digits) (+ (* i digit-size) j))
+ (declare (type (mod #.sb!vm:n-word-bits) j))))))))
+
+;;; Multiply a bignum buffer with a fixnum or a digit, storing the
+;;; result in another bignum buffer, and without using any
+;;; temporaries. Inlined to avoid boxing smallnum if it's actually a
+;;; digit. Needed by GCD, should possibly OAOO with
+;;; MULTIPLY-BIGNUM-AND-FIXNUM.
+(declaim (inline multiply-bignum-buffer-and-smallnum-to-buffer))
+(defun multiply-bignum-buffer-and-smallnum-to-buffer (bignum bignum-len
+ smallnum res)
+ (declare (type bignum-type bignum))
+ (let* ((bignum-plus-p (%bignum-0-or-plusp bignum bignum-len))
+ (smallnum-plus-p (not (minusp smallnum)))
+ (smallnum (if smallnum-plus-p smallnum (- smallnum)))
+ (carry-digit 0))
+ (declare (type bignum-type bignum res)
+ (type bignum-index bignum-len)
+ (type bignum-element-type smallnum carry-digit))
+ (unless bignum-plus-p
+ (negate-bignum-buffer-in-place bignum bignum-len))
+ (dotimes (index bignum-len)
+ (declare (type bignum-index index))
+ (multiple-value-bind (next-digit low)
+ (%multiply-and-add (%bignum-ref bignum index)
+ smallnum
+ carry-digit)
+ (declare (type bignum-element-type next-digit low))
+ (setf carry-digit next-digit)
+ (setf (%bignum-ref res index) low)))
+ (setf (%bignum-ref res bignum-len) carry-digit)
+ (unless bignum-plus-p
+ (negate-bignum-buffer-in-place bignum bignum-len))
+ (let ((res-len (%normalize-bignum-buffer res (1+ bignum-len))))
+ (unless (eq bignum-plus-p smallnum-plus-p)
+ (negate-bignum-buffer-in-place res res-len))
+ res-len)))
+
+;;; Given U and V, return U / V mod 2^32. Implements the algorithm in the
+;;; paper, but uses some clever bit-twiddling nicked from Nickle to do it.
+(declaim (inline bmod))
+(defun bmod (u v)
+ (let ((ud (%bignum-ref u 0))
+ (vd (%bignum-ref v 0))
+ (umask 0)
+ (imask 1)
+ (m 0))
+ (declare (type (unsigned-byte #.sb!vm:n-word-bits) ud vd umask imask m))
+ (dotimes (i digit-size)
+ (setf umask (logior umask imask))
+ (unless (zerop (logand ud umask))
+ (setf ud (M (- ud vd)))
+ (setf m (M (logior m imask))))
+ (setf imask (M (ash imask 1)))
+ (setf vd (M (ash vd 1))))
+ m))
+
+(defun dmod (u u-len v v-len tmp1)
+ (loop while (> (bignum-buffer-integer-length u u-len)
+ (+ (bignum-buffer-integer-length v v-len)
+ digit-size))
+ do
+ (unless (zerop (%bignum-ref u 0))
+ (let* ((bmod (bmod u v))
+ (tmp1-len (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
+ bmod
+ tmp1)))
+ (setf u-len (subtract-bignum-buffers u u-len
+ tmp1 tmp1-len
+ u))
+ (bignum-abs-buffer u u-len)))
+ (gcd-assert (zerop (%bignum-ref u 0)))
+ (setf u-len (bignum-buffer-ashift-right u u-len digit-size)))
+ (let* ((d (+ 1 (- (bignum-buffer-integer-length u u-len)
+ (bignum-buffer-integer-length v v-len))))
+ (n (1- (ash 1 d))))
+ (declare (type (unsigned-byte #.(integer-length #.sb!vm:n-word-bits)) d)
+ (type (unsigned-byte #.sb!vm:n-word-bits) n))
+ (gcd-assert (>= d 0))
+ (unless (zerop (logand (%bignum-ref u 0) n))
+ (let ((tmp1-len
+ (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
+ (logand n (bmod u
+ v))
+ tmp1)))
+ (setf u-len (subtract-bignum-buffers u u-len
+ tmp1 tmp1-len
+ u))
+ (bignum-abs-buffer u u-len)))
+ u-len))
+
+(defconstant lower-ones-digit (1- (ash 1 (truncate sb!vm:n-word-bits 2))))
+
+;;; Find D and N such that (LOGAND ALL-ONES-DIGIT (- (* D X) (* N Y))) is 0,
+;;; (< 0 N LOWER-ONES-DIGIT) and (< 0 (ABS D) LOWER-ONES-DIGIT).
+(defun reduced-ratio-mod (x y)
+ (let* ((c (bmod x y))
+ (n1 c)
+ (d1 1)
+ (n2 (M (1+ (M (lognot n1)))))
+ (d2 (M -1)))
+ (declare (type (unsigned-byte #.sb!vm:n-word-bits) n1 d1 n2 d2))
+ (loop while (> n2 (expt 2 (truncate digit-size 2))) do
+ (loop for i of-type (mod #.sb!vm:n-word-bits)
+ downfrom (- (integer-length n1) (integer-length n2))
+ while (>= n1 n2) do
+ (when (>= n1 (M (ash n2 i)))
+ (psetf n1 (M (- n1 (M (ash n2 i))))
+ d1 (M (- d1 (M (ash d2 i)))))))
+ (psetf n1 n2
+ d1 d2
+ n2 n1
+ d2 d1))
+ (values n2 (if (>= d2 (expt 2 (1- digit-size)))
+ (lognot (logand most-positive-fixnum (lognot d2)))
+ (logand lower-ones-digit d2)))))
+
+
+(defun copy-bignum (a &optional (len (%bignum-length a)))
+ (let ((b (%allocate-bignum len)))
+ (bignum-replace b a)
+ (%bignum-set-length b len)
+ b))
+
+;;; Allocate a single word bignum that holds fixnum. This is useful when
+;;; we are trying to mix fixnum and bignum operands.
+#!-sb-fluid (declaim (inline make-small-bignum))
+(defun make-small-bignum (fixnum)
+ (let ((res (%allocate-bignum 1)))
+ (setf (%bignum-ref res 0) (%fixnum-to-digit fixnum))
+ res))
+
+;; When the larger number is less than this many bignum digits long, revert
+;; to old algorithm.
+(defparameter *accelerated-gcd-cutoff* 3)
+
+;;; Alternate between k-ary reduction with the help of
+;;; REDUCED-RATIO-MOD and digit modulus reduction via DMOD. Once the
+;;; arguments get small enough, drop through to BIGNUM-MOD-GCD (since
+;;; k-ary reduction can introduce spurious factors, which need to be
+;;; filtered out). Reference: Kenneth Weber, "The accelerated integer
+;;; GCD algorithm", ACM Transactions on Mathematical Software, volume
+;;; 21, number 1, March 1995, epp. 111-122.
+(defun bignum-gcd (u0 v0)
+ (declare (type bignum-type u0 v0))
+ (let* ((u1 (if (%bignum-0-or-plusp u0 (%bignum-length u0))
+ u0
+ (negate-bignum u0 nil)))
+ (v1 (if (%bignum-0-or-plusp v0 (%bignum-length v0))
+ v0
+ (negate-bignum v0 nil))))
+ (if (zerop v1)
+ (return-from bignum-gcd u1))
+ (when (> u1 v1)
+ (rotatef u1 v1))
+ (let ((n (mod v1 u1)))
+ (setf v1 (if (fixnump n)
+ (make-small-bignum n)
+ n)))
+ (if (and (= 1 (%bignum-length v1))
+ (zerop (%bignum-ref v1 0)))
+ (return-from bignum-gcd (%normalize-bignum u1
+ (%bignum-length u1))))
+ (let* ((buffer-len (+ 2 (%bignum-length u1)))
+ (u (%allocate-bignum buffer-len))
+ (u-len (%bignum-length u1))
+ (v (%allocate-bignum buffer-len))
+ (v-len (%bignum-length v1))
+ (tmp1 (%allocate-bignum buffer-len))
+ (tmp1-len 0)
+ (tmp2 (%allocate-bignum buffer-len))
+ (tmp2-len 0)
+ (factors-of-two
+ (bignum-factors-of-two u1 (%bignum-length u1)
+ v1 (%bignum-length v1))))
+ (declare (type (or null bignum-index)
+ buffer-len u-len v-len tmp1-len tmp2-len))
+ (bignum-replace u u1)
+ (bignum-replace v v1)
+ (setf u-len
+ (make-gcd-bignum-odd u
+ (bignum-buffer-ashift-right u u-len
+ factors-of-two)))
+ (setf v-len
+ (make-gcd-bignum-odd v
+ (bignum-buffer-ashift-right v v-len
+ factors-of-two)))
+ (loop until (or (< u-len *accelerated-gcd-cutoff*)
+ (not v-len)
+ (zerop v-len)
+ (and (= 1 v-len)
+ (zerop (%bignum-ref v 0))))
+ do
+ (gcd-assert (= buffer-len (%bignum-length u)
+ (%bignum-length v)
+ (%bignum-length tmp1)
+ (%bignum-length tmp2)))
+ (if (> (bignum-buffer-integer-length u u-len)
+ (+ #.(truncate sb!vm:n-word-bits 4)
+ (bignum-buffer-integer-length v v-len)))
+ (setf u-len (dmod u u-len
+ v v-len
+ tmp1))
+ (multiple-value-bind (n d) (reduced-ratio-mod u v)
+ (setf tmp1-len
+ (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
+ n tmp1))
+ (setf tmp2-len
+ (multiply-bignum-buffer-and-smallnum-to-buffer u u-len
+ d tmp2))
+ (gcd-assert (= (copy-bignum tmp2 tmp2-len)
+ (* (copy-bignum u u-len) d)))
+ (gcd-assert (= (copy-bignum tmp1 tmp1-len)
+ (* (copy-bignum v v-len) n)))
+ (setf u-len
+ (subtract-bignum-buffers-with-len tmp1 tmp1-len
+ tmp2 tmp2-len
+ u
+ (1+ (max tmp1-len
+ tmp2-len))))
+ (gcd-assert (or (zerop (- (copy-bignum tmp1 tmp1-len)
+ (copy-bignum tmp2 tmp2-len)))
+ (= (copy-bignum u u-len)
+ (- (copy-bignum tmp1 tmp1-len)
+ (copy-bignum tmp2 tmp2-len)))))
+ (bignum-abs-buffer u u-len)
+ (gcd-assert (zerop (M u)))))
+ (setf u-len (make-gcd-bignum-odd u u-len))
+ (rotatef u v)
+ (rotatef u-len v-len))
+ (setf u (copy-bignum u u-len))
+ (let ((n (bignum-mod-gcd v1 u)))
+ (ash (bignum-mod-gcd u1 (if (fixnump n)
+ (make-small-bignum n)
+ n))
+ factors-of-two)))))
+
+(defun bignum-mod-gcd (a b)
(declare (type bignum-type a b))
- (let* ((a (if (%bignum-0-or-plusp a (%bignum-length a))
- a
- (negate-bignum a nil)))
- (b (if (%bignum-0-or-plusp b (%bignum-length b))
- b
- (negate-bignum b nil)))
- (len-a (%bignum-length a))
+ (when (< a b)
+ (rotatef a b))
+ ;; While the length difference of A and B is sufficiently large,
+ ;; reduce using MOD (slowish, but it should equalize the sizes of
+ ;; A and B pretty quickly). After that, use the binary GCD
+ ;; algorithm to handle the rest.
+ (loop until (and (= (%bignum-length b) 1) (zerop (%bignum-ref b 0))) do
+ (when (<= (%bignum-length a) (1+ (%bignum-length b)))
+ (return-from bignum-mod-gcd (bignum-binary-gcd a b)))
+ (let ((rem (mod a b)))
+ (if (fixnump rem)
+ (setf a (make-small-bignum rem))
+ (setf a rem))
+ (rotatef a b)))
+ (if (= (%bignum-length a) 1)
+ (%normalize-bignum a 1)
+ a))
+
+(defun bignum-binary-gcd (a b)
+ (declare (type bignum-type a b))
+ (let* ((len-a (%bignum-length a))
(len-b (%bignum-length b)))
- (declare (type bignum-index len-a len-b))
+ (declare (type bignum-index len-a len-b))
(with-bignum-buffers ((a-buffer len-a a)
(b-buffer len-b b)
(res-buffer (max len-a len-b)))
((%digit-greater a-digit b-digit)
(return
(values a b len-b res
- (subtract-bignum-buffers a len-a b len-b res))))
+ (subtract-bignum-buffers a len-a b len-b
+ res))))
(t
(return
(values b a len-a res
(do ((digit (%bignum-ref a index) (%ashr digit 1))
(increment 0 (1+ increment)))
((zerop digit))
- (declare (type (mod 32) increment))
+ (declare (type (mod #.sb!vm:n-word-bits) increment))
(when (oddp digit)
(return-from make-gcd-bignum-odd
(bignum-buffer-ashift-right a len-a
(+ (* index digit-size)
increment)))))))
-(defun bignum-factors-of-two (a len-a b len-b)
- (declare (type bignum-index len-a len-b) (type bignum-type a))
- (do ((i 0 (1+ i))
- (end (min len-a len-b)))
- ((= i end) (error "Unexpected zero bignums?"))
- (declare (type bignum-index i end))
- (let ((or-digits (%logior (%bignum-ref a i) (%bignum-ref b i))))
- (unless (zerop or-digits)
- (return (do ((j 0 (1+ j))
- (or-digits or-digits (%ashr or-digits 1)))
- ((oddp or-digits) (+ (* i digit-size) j))
- (declare (type (mod 32) j))))))))
\f
;;;; negation
;;; This assumes bignum is positive; that is, the result of negating it will
;;; stay in the provided allocated bignum.
-(defun negate-bignum-in-place (bignum)
- (bignum-negate-loop bignum (%bignum-length bignum) bignum)
+(defun negate-bignum-buffer-in-place (bignum bignum-len)
+ (bignum-negate-loop bignum bignum-len bignum)
bignum)
+
+(defun negate-bignum-in-place (bignum)
+ (declare (inline negate-bignum-buffer-in-place))
+ (negate-bignum-buffer-in-place bignum (%bignum-length bignum)))
+
+(defun bignum-abs-buffer (bignum len)
+ (unless (%bignum-0-or-plusp bignum len)
+ (negate-bignum-in-place bignum len)))
\f
;;;; shifting
-(defconstant all-ones-digit #xFFFFFFFF)
-
(eval-when (:compile-toplevel :execute)
;;; This macro is used by BIGNUM-ASHIFT-RIGHT, BIGNUM-BUFFER-ASHIFT-RIGHT, and
(%normalize-bignum res res-len))
res)))))
((> count bignum-len)
- 0)
+ (if (%bignum-0-or-plusp bignum bignum-len) 0 -1))
;; Since a FIXNUM should be big enough to address anything in
;; memory, including arrays of bits, and since arrays of bits
;; take up about the same space as corresponding fixnums, there
;; should be no way that we fall through to this case: any shift
;; right by a bignum should give zero. But let's check anyway:
- (t (error "bignum overflow: can't shift right by ~S")))))
+ (t (error "bignum overflow: can't shift right by ~S" count)))))
(defun bignum-ashift-right-digits (bignum digits)
(declare (type bignum-type bignum)
(bignum-ashift-left-unaligned bignum digits n-bits res-len))))
;; Left shift by a number too big to be represented as a fixnum
;; would exceed our memory capacity, since a fixnum is big enough
- ;; index any array, including a bit array.
+ ;; to index any array, including a bit array.
(error "can't represent result of left shift")))
(defun bignum-ashift-left-digits (bignum bignum-len digits)
(declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
(let ((res (dpb exp
sb!vm:single-float-exponent-byte
- (logandc2 (sb!ext:truly-the (unsigned-byte 31)
+ (logandc2 (sb!ext:truly-the (unsigned-byte #.(1- sb!vm:n-word-bits))
(%bignum-ref bits 1))
sb!vm:single-float-hidden-bit))))
(make-single-float
(declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
(let ((hi (dpb exp
sb!vm:double-float-exponent-byte
- (logandc2 (sb!ext:truly-the (unsigned-byte 31)
+ (logandc2 (sb!ext:truly-the (unsigned-byte #.(1- sb!vm:n-word-bits))
(%bignum-ref bits 2))
sb!vm:double-float-hidden-bit))))
(make-double-float
(declare (type bignum-index len))
(let ((exp (+ exp bias)))
(when (> exp max)
- (error "Too large to be represented as a ~S:~% ~S"
- format x))
+ ;; Why a SIMPLE-TYPE-ERROR? Well, this is mainly
+ ;; called by COERCE, which requires an error of
+ ;; TYPE-ERROR if the conversion can't happen
+ ;; (except in certain circumstances when we are
+ ;; coercing to a FUNCTION) -- CSR, 2002-09-18
+ (error 'simple-type-error
+ :format-control "Too large to be represented as a ~S:~% ~S"
+ :format-arguments (list format x)
+ :expected-type format
+ :datum x))
exp)))
(cond
(t
(round-up))))))
\f
-;;;; integer length and logcount
+;;;; integer length and logbitp/logcount
-(defun bignum-integer-length (bignum)
+(defun bignum-buffer-integer-length (bignum len)
(declare (type bignum-type bignum))
- (let* ((len (%bignum-length bignum))
- (len-1 (1- len))
+ (let* ((len-1 (1- len))
(digit (%bignum-ref bignum len-1)))
(declare (type bignum-index len len-1)
(type bignum-element-type digit))
(+ (integer-length (%fixnum-digit-with-correct-sign digit))
(* len-1 digit-size))))
+(defun bignum-integer-length (bignum)
+ (declare (type bignum-type bignum))
+ (bignum-buffer-integer-length bignum (%bignum-length bignum)))
+
+(defun bignum-logbitp (index bignum)
+ (declare (type bignum-type bignum))
+ (let ((len (%bignum-length bignum)))
+ (declare (type bignum-index len))
+ (multiple-value-bind (word-index bit-index)
+ (floor index digit-size)
+ (if (>= word-index len)
+ (not (bignum-plus-p bignum))
+ (not (zerop (logand (%bignum-ref bignum word-index)
+ (ash 1 bit-index))))))))
+
(defun bignum-logcount (bignum)
(declare (type bignum-type bignum))
(let* ((length (%bignum-length bignum))
FOR NOW WE DON'T USE LDB OR DPB. WE USE SHIFTS AND MASKS IN NUMBERS.LISP WHICH
IS LESS EFFICIENT BUT EASIER TO MAINTAIN. BILL SAYS THIS CODE CERTAINLY WORKS!
-(defconstant maximum-fixnum-bits #!+ibm-rt-pc 27 #!-ibm-rt-pc 30)
+(defconstant maximum-fixnum-bits (- sb!vm:n-word-bits sb!vm:n-lowtag-bits))
(defun bignum-load-byte (byte bignum)
(declare (type bignum-type bignum))
(- pos))
;; LOGAND should be unnecessary here
;; with a logical right shift or a
- ;; correct unsigned-byte-32 one.
+ ;; correct digit-sized one.
(%make-ones available-bits))))
(if (%bignum-0-or-plusp bignum bignum-len)
res
(- pos))
;; LOGAND should be unnecessary here with
;; a logical right shift or a correct
- ;; unsigned-byte-32 one.
+ ;; digit-sized one.
high-mask))))))))))
;;; This returns a bignum result of loading a byte from a bignum. In order, we
(setf (%bignum-ref res j)
(logand (ash (%bignum-ref bignum i) minus-start-pos)
;; LOGAND should be unnecessary here with a logical
- ;; right shift or a correct unsigned-byte-32 one.
+ ;; right shift or a correct digit-sized one.
high-mask))
(when (%bignum-0-or-plusp bignum bignum-len)
(setf (%bignum-ref res j)
(setf (%bignum-ref res j)
(logand (ash (%bignum-ref bignum i) minus-start-pos)
;; LOGAND should be unnecessary here with a logical
- ;; right shift or a correct unsigned-byte-32 one.
+ ;; right shift or a correct digit-sized one.
high-mask))
(unless (%bignum-0-or-plusp bignum bignum-len)
;; Fill in upper half of this result digit with 1's.
(setf (%bignum-ref result j+1)
(logand (ash digit minus-high-bits)
;; LOGAND should be unnecessary here with a logical right
- ;; shift or a correct unsigned-byte-32 one.
+ ;; shift or a correct digit-sized one.
low-mask))))))
|#
\f
;;; j = last digit of y.
;;;
;;; compute guess.
-;;; if x[i] = y[j] then g = #xFFFFFFFF
+;;; if x[i] = y[j] then g = (1- (ash 1 digit-size))
;;; else g = x[i]x[i-1]/y[j].
;;;
;;; check guess.
;;; %UNSIGNED-MULTIPLY returns b and c defined below.
;;; a = x[i-1] - (logand (* g y[j]) #xFFFFFFFF).
;;; Use %UNSIGNED-MULTIPLY taking low-order result.
-;;; b = (logand (ash (* g y[j-1]) -32) #xFFFFFFFF).
-;;; c = (logand (* g y[j-1]) #xFFFFFFFF).
+;;; b = (logand (ash (* g y[j-1]) (- digit-size)) (1- (ash 1 digit-size))).
+;;; c = (logand (* g y[j-1]) (1- (ash 1 digit-size))).
;;; if a < b, okay.
;;; if a > b, guess is too high
;;; g = g - 1; go back to "check guess".
;;;
;;; Normalize quotient and remainder. Cons result if necessary.
-;;; These are used by BIGNUM-TRUNCATE and friends in the general case.
-(defvar *truncate-x*)
-(defvar *truncate-y*)
-;;; This divides x by y returning the quotient and remainder. In the general
-;;; case, we shift y to setup for the algorithm, and we use two buffers to save
-;;; consing intermediate values. X gets destructively modified to become the
-;;; remainder, and we have to shift it to account for the initial Y shift.
-;;; After we multiple bind q and r, we first fix up the signs and then return
-;;; the normalized results.
+;;; This used to be split into multiple functions, which shared state
+;;; in special variables *TRUNCATE-X* and *TRUNCATE-Y*. Having so many
+;;; special variable accesses in tight inner loops was having a large
+;;; effect on performance, so the helper functions have now been
+;;; refactored into local functions and the special variables into
+;;; lexicals. There was also a lot of boxing and unboxing of
+;;; (UNSIGNED-BYTE 32)'s going on, which this refactoring
+;;; eliminated. This improves the performance on some CL-BENCH tests
+;;; by up to 50%, which is probably signigicant enough to justify the
+;;; reduction in readability that was introduced. --JES, 2004-08-07
(defun bignum-truncate (x y)
(declare (type bignum-type x y))
- (let* ((x-plusp (%bignum-0-or-plusp x (%bignum-length x)))
- (y-plusp (%bignum-0-or-plusp y (%bignum-length y)))
- (x (if x-plusp x (negate-bignum x nil)))
- (y (if y-plusp y (negate-bignum y nil)))
- (len-x (%bignum-length x))
- (len-y (%bignum-length y)))
- (multiple-value-bind (q r)
- (cond ((< len-y 2)
- (bignum-truncate-single-digit x len-x y))
- ((plusp (bignum-compare y x))
- (let ((res (%allocate-bignum len-x)))
- (dotimes (i len-x)
- (setf (%bignum-ref res i) (%bignum-ref x i)))
- (values 0 res)))
- (t
- (let ((len-x+1 (1+ len-x)))
- (with-bignum-buffers ((*truncate-x* len-x+1)
- (*truncate-y* (1+ len-y)))
- (let ((y-shift (shift-y-for-truncate y)))
- (shift-and-store-truncate-buffers x len-x y len-y y-shift)
- (values (do-truncate len-x+1 len-y)
- ;; DO-TRUNCATE must execute first.
- (cond
- ((zerop y-shift)
- (let ((res (%allocate-bignum len-y)))
- (declare (type bignum-type res))
- (bignum-replace res *truncate-x* :end2 len-y)
- (%normalize-bignum res len-y)))
- (t
- (shift-right-unaligned
- *truncate-x* 0 y-shift len-y
- ((= j res-len-1)
- (setf (%bignum-ref res j)
- (%ashr (%bignum-ref *truncate-x* i)
- y-shift))
- (%normalize-bignum res res-len))
- res)))))))))
- (let ((quotient (cond ((eq x-plusp y-plusp) q)
- ((typep q 'fixnum) (the fixnum (- q)))
- (t (negate-bignum-in-place q))))
- (rem (cond (x-plusp r)
- ((typep r 'fixnum) (the fixnum (- r)))
- (t (negate-bignum-in-place r)))))
- (values (if (typep quotient 'fixnum)
- quotient
- (%normalize-bignum quotient (%bignum-length quotient)))
- (if (typep rem 'fixnum)
- rem
- (%normalize-bignum rem (%bignum-length rem))))))))
-
-;;; This divides x by y when y is a single bignum digit. BIGNUM-TRUNCATE fixes
-;;; up the quotient and remainder with respect to sign and normalization.
-;;;
-;;; We don't have to worry about shifting y to make its most significant digit
-;;; sufficiently large for %FLOOR to return 32-bit quantities for the q-digit
-;;; and r-digit. If y is a single digit bignum, it is already large enough
-;;; for %FLOOR. That is, it has some bits on pretty high in the digit.
-(defun bignum-truncate-single-digit (x len-x y)
- (declare (type bignum-index len-x))
- (let ((q (%allocate-bignum len-x))
- (r 0)
- (y (%bignum-ref y 0)))
- (declare (type bignum-element-type r y))
- (do ((i (1- len-x) (1- i)))
- ((minusp i))
- (multiple-value-bind (q-digit r-digit) (%floor r (%bignum-ref x i) y)
- (declare (type bignum-element-type q-digit r-digit))
- (setf (%bignum-ref q i) q-digit)
- (setf r r-digit)))
- (let ((rem (%allocate-bignum 1)))
- (setf (%bignum-ref rem 0) r)
- (values q rem))))
-
-;;; This divides *truncate-x* by *truncate-y*, and len-x and len-y tell us how
-;;; much of the buffers we care about. TRY-BIGNUM-TRUNCATE-GUESS modifies
-;;; *truncate-x* on each interation, and this buffer becomes our remainder.
-;;;
-;;; *truncate-x* definitely has at least three digits, and it has one more than
-;;; *truncate-y*. This keeps i, i-1, i-2, and low-x-digit happy. Thanks to
-;;; SHIFT-AND-STORE-TRUNCATE-BUFFERS.
-(defun do-truncate (len-x len-y)
- (declare (type bignum-index len-x len-y))
- (let* ((len-q (- len-x len-y))
- ;; Add one for extra sign digit in case high bit is on.
- (q (%allocate-bignum (1+ len-q)))
- (k (1- len-q))
- (y1 (%bignum-ref *truncate-y* (1- len-y)))
- (y2 (%bignum-ref *truncate-y* (- len-y 2)))
- (i (1- len-x))
- (i-1 (1- i))
- (i-2 (1- i-1))
- (low-x-digit (- i len-y)))
- (declare (type bignum-index len-q k i i-1 i-2 low-x-digit)
- (type bignum-element-type y1 y2))
- (loop
- (setf (%bignum-ref q k)
- (try-bignum-truncate-guess
- ;; This modifies *truncate-x*. Must access elements each pass.
- (bignum-truncate-guess y1 y2
- (%bignum-ref *truncate-x* i)
- (%bignum-ref *truncate-x* i-1)
- (%bignum-ref *truncate-x* i-2))
- len-y low-x-digit))
- (cond ((zerop k) (return))
- (t (decf k)
- (decf low-x-digit)
- (shiftf i i-1 i-2 (1- i-2)))))
- q))
-
-;;; This takes a digit guess, multiplies it by *truncate-y* for a result one
-;;; greater in length than len-y, and subtracts this result from *truncate-x*.
-;;; Low-x-digit is the first digit of x to start the subtraction, and we know x
-;;; is long enough to subtract a len-y plus one length bignum from it. Next we
-;;; check the result of the subtraction, and if the high digit in x became
-;;; negative, then our guess was one too big. In this case, return one less
-;;; than guess passed in, and add one value of y back into x to account for
-;;; subtracting one too many. Knuth shows that the guess is wrong on the order
-;;; of 3/b, where b is the base (2 to the digit-size power) -- pretty rarely.
-(defun try-bignum-truncate-guess (guess len-y low-x-digit)
- (declare (type bignum-index low-x-digit len-y)
- (type bignum-element-type guess))
- (let ((carry-digit 0)
- (borrow 1)
- (i low-x-digit))
- (declare (type bignum-element-type carry-digit)
- (type bignum-index i)
- (fixnum borrow))
- ;; Multiply guess and divisor, subtracting from dividend simultaneously.
- (dotimes (j len-y)
- (multiple-value-bind (high-digit low-digit)
- (%multiply-and-add guess
- (%bignum-ref *truncate-y* j)
- carry-digit)
- (declare (type bignum-element-type high-digit low-digit))
- (setf carry-digit high-digit)
- (multiple-value-bind (x temp-borrow)
- (%subtract-with-borrow (%bignum-ref *truncate-x* i)
- low-digit
- borrow)
- (declare (type bignum-element-type x)
- (fixnum temp-borrow))
- (setf (%bignum-ref *truncate-x* i) x)
- (setf borrow temp-borrow)))
- (incf i))
- (setf (%bignum-ref *truncate-x* i)
- (%subtract-with-borrow (%bignum-ref *truncate-x* i)
- carry-digit borrow))
- ;; See whether guess is off by one, adding one Y back in if necessary.
- (cond ((%digit-0-or-plusp (%bignum-ref *truncate-x* i))
- guess)
- (t
- ;; If subtraction has negative result, add one divisor value back
- ;; in. The guess was one too large in magnitude.
- (let ((i low-x-digit)
- (carry 0))
+ (let (truncate-x truncate-y)
+ (labels
+ ;;; Divide X by Y when Y is a single bignum digit. BIGNUM-TRUNCATE
+ ;;; fixes up the quotient and remainder with respect to sign and
+ ;;; normalization.
+ ;;;
+ ;;; We don't have to worry about shifting Y to make its most
+ ;;; significant digit sufficiently large for %FLOOR to return
+ ;;; digit-size quantities for the q-digit and r-digit. If Y is
+ ;;; a single digit bignum, it is already large enough for
+ ;;; %FLOOR. That is, it has some bits on pretty high in the
+ ;;; digit.
+ ((bignum-truncate-single-digit (x len-x y)
+ (declare (type bignum-index len-x))
+ (let ((q (%allocate-bignum len-x))
+ (r 0)
+ (y (%bignum-ref y 0)))
+ (declare (type bignum-element-type r y))
+ (do ((i (1- len-x) (1- i)))
+ ((minusp i))
+ (multiple-value-bind (q-digit r-digit)
+ (%floor r (%bignum-ref x i) y)
+ (declare (type bignum-element-type q-digit r-digit))
+ (setf (%bignum-ref q i) q-digit)
+ (setf r r-digit)))
+ (let ((rem (%allocate-bignum 1)))
+ (setf (%bignum-ref rem 0) r)
+ (values q rem))))
+ ;;; This returns a guess for the next division step. Y1 is the
+ ;;; highest y digit, and y2 is the second to highest y
+ ;;; digit. The x... variables are the three highest x digits
+ ;;; for the next division step.
+ ;;;
+ ;;; From Knuth, our guess is either all ones or x-i and x-i-1
+ ;;; divided by y1, depending on whether x-i and y1 are the
+ ;;; same. We test this guess by determining whether guess*y2
+ ;;; is greater than the three high digits of x minus guess*y1
+ ;;; shifted left one digit:
+ ;;; ------------------------------
+ ;;; | x-i | x-i-1 | x-i-2 |
+ ;;; ------------------------------
+ ;;; ------------------------------
+ ;;; - | g*y1 high | g*y1 low | 0 |
+ ;;; ------------------------------
+ ;;; ... < guess*y2 ???
+ ;;; If guess*y2 is greater, then we decrement our guess by one
+ ;;; and try again. This returns a guess that is either
+ ;;; correct or one too large.
+ (bignum-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
+ (declare (type bignum-element-type y1 y2 x-i x-i-1 x-i-2))
+ (let ((guess (if (%digit-compare x-i y1)
+ all-ones-digit
+ (%floor x-i x-i-1 y1))))
+ (declare (type bignum-element-type guess))
+ (loop
+ (multiple-value-bind (high-guess*y1 low-guess*y1)
+ (%multiply guess y1)
+ (declare (type bignum-element-type low-guess*y1
+ high-guess*y1))
+ (multiple-value-bind (high-guess*y2 low-guess*y2)
+ (%multiply guess y2)
+ (declare (type bignum-element-type high-guess*y2
+ low-guess*y2))
+ (multiple-value-bind (middle-digit borrow)
+ (%subtract-with-borrow x-i-1 low-guess*y1 1)
+ (declare (type bignum-element-type middle-digit)
+ (fixnum borrow))
+ ;; Supplying borrow of 1 means there was no
+ ;; borrow, and we know x-i-2 minus 0 requires
+ ;; no borrow.
+ (let ((high-digit (%subtract-with-borrow x-i
+ high-guess*y1
+ borrow)))
+ (declare (type bignum-element-type high-digit))
+ (if (and (%digit-compare high-digit 0)
+ (or (%digit-greater high-guess*y2
+ middle-digit)
+ (and (%digit-compare middle-digit
+ high-guess*y2)
+ (%digit-greater low-guess*y2
+ x-i-2))))
+ (setf guess (%subtract-with-borrow guess 1 1))
+ (return guess)))))))))
+ ;;; Divide TRUNCATE-X by TRUNCATE-Y, returning the quotient
+ ;;; and destructively modifying TRUNCATE-X so that it holds
+ ;;; the remainder.
+ ;;;
+ ;;; LEN-X and LEN-Y tell us how much of the buffers we care about.
+ ;;;
+ ;;; TRUNCATE-X definitely has at least three digits, and it has one
+ ;;; more than TRUNCATE-Y. This keeps i, i-1, i-2, and low-x-digit
+ ;;; happy. Thanks to SHIFT-AND-STORE-TRUNCATE-BUFFERS.
+ (return-quotient-leaving-remainder (len-x len-y)
+ (declare (type bignum-index len-x len-y))
+ (let* ((len-q (- len-x len-y))
+ ;; Add one for extra sign digit in case high bit is on.
+ (q (%allocate-bignum (1+ len-q)))
+ (k (1- len-q))
+ (y1 (%bignum-ref truncate-y (1- len-y)))
+ (y2 (%bignum-ref truncate-y (- len-y 2)))
+ (i (1- len-x))
+ (i-1 (1- i))
+ (i-2 (1- i-1))
+ (low-x-digit (- i len-y)))
+ (declare (type bignum-index len-q k i i-1 i-2 low-x-digit)
+ (type bignum-element-type y1 y2))
+ (loop
+ (setf (%bignum-ref q k)
+ (try-bignum-truncate-guess
+ ;; This modifies TRUNCATE-X. Must access
+ ;; elements each pass.
+ (bignum-truncate-guess y1 y2
+ (%bignum-ref truncate-x i)
+ (%bignum-ref truncate-x i-1)
+ (%bignum-ref truncate-x i-2))
+ len-y low-x-digit))
+ (cond ((zerop k) (return))
+ (t (decf k)
+ (decf low-x-digit)
+ (shiftf i i-1 i-2 (1- i-2)))))
+ q))
+ ;;; This takes a digit guess, multiplies it by TRUNCATE-Y for a
+ ;;; result one greater in length than LEN-Y, and subtracts this result
+ ;;; from TRUNCATE-X. LOW-X-DIGIT is the first digit of X to start
+ ;;; the subtraction, and we know X is long enough to subtract a LEN-Y
+ ;;; plus one length bignum from it. Next we check the result of the
+ ;;; subtraction, and if the high digit in X became negative, then our
+ ;;; guess was one too big. In this case, return one less than GUESS
+ ;;; passed in, and add one value of Y back into X to account for
+ ;;; subtracting one too many. Knuth shows that the guess is wrong on
+ ;;; the order of 3/b, where b is the base (2 to the digit-size power)
+ ;;; -- pretty rarely.
+ (try-bignum-truncate-guess (guess len-y low-x-digit)
+ (declare (type bignum-index low-x-digit len-y)
+ (type bignum-element-type guess))
+ (let ((carry-digit 0)
+ (borrow 1)
+ (i low-x-digit))
+ (declare (type bignum-element-type carry-digit)
+ (type bignum-index i)
+ (fixnum borrow))
+ ;; Multiply guess and divisor, subtracting from dividend
+ ;; simultaneously.
(dotimes (j len-y)
- (multiple-value-bind (v k)
- (%add-with-carry (%bignum-ref *truncate-y* j)
- (%bignum-ref *truncate-x* i)
- carry)
- (declare (type bignum-element-type v))
- (setf (%bignum-ref *truncate-x* i) v)
- (setf carry k))
+ (multiple-value-bind (high-digit low-digit)
+ (%multiply-and-add guess
+ (%bignum-ref truncate-y j)
+ carry-digit)
+ (declare (type bignum-element-type high-digit low-digit))
+ (setf carry-digit high-digit)
+ (multiple-value-bind (x temp-borrow)
+ (%subtract-with-borrow (%bignum-ref truncate-x i)
+ low-digit
+ borrow)
+ (declare (type bignum-element-type x)
+ (fixnum temp-borrow))
+ (setf (%bignum-ref truncate-x i) x)
+ (setf borrow temp-borrow)))
(incf i))
- (setf (%bignum-ref *truncate-x* i)
- (%add-with-carry (%bignum-ref *truncate-x* i) 0 carry)))
- (%subtract-with-borrow guess 1 1)))))
+ (setf (%bignum-ref truncate-x i)
+ (%subtract-with-borrow (%bignum-ref truncate-x i)
+ carry-digit borrow))
+ ;; See whether guess is off by one, adding one
+ ;; Y back in if necessary.
+ (cond ((%digit-0-or-plusp (%bignum-ref truncate-x i))
+ guess)
+ (t
+ ;; If subtraction has negative result, add one
+ ;; divisor value back in. The guess was one too
+ ;; large in magnitude.
+ (let ((i low-x-digit)
+ (carry 0))
+ (dotimes (j len-y)
+ (multiple-value-bind (v k)
+ (%add-with-carry (%bignum-ref truncate-y j)
+ (%bignum-ref truncate-x i)
+ carry)
+ (declare (type bignum-element-type v))
+ (setf (%bignum-ref truncate-x i) v)
+ (setf carry k))
+ (incf i))
+ (setf (%bignum-ref truncate-x i)
+ (%add-with-carry (%bignum-ref truncate-x i)
+ 0 carry)))
+ (%subtract-with-borrow guess 1 1)))))
+ ;;; This returns the amount to shift y to place a one in the
+ ;;; second highest bit. Y must be positive. If the last digit
+ ;;; of y is zero, then y has a one in the previous digit's
+ ;;; sign bit, so we know it will take one less than digit-size
+ ;;; to get a one where we want. Otherwise, we count how many
+ ;;; right shifts it takes to get zero; subtracting this value
+ ;;; from digit-size tells us how many high zeros there are
+ ;;; which is one more than the shift amount sought.
+ ;;;
+ ;;; Note: This is exactly the same as one less than the
+ ;;; integer-length of the last digit subtracted from the
+ ;;; digit-size.
+ ;;;
+ ;;; We shift y to make it sufficiently large that doing the
+ ;;; 2*digit-size by digit-size %FLOOR calls ensures the quotient and
+ ;;; remainder fit in digit-size.
+ (shift-y-for-truncate (y)
+ (let* ((len (%bignum-length y))
+ (last (%bignum-ref y (1- len))))
+ (declare (type bignum-index len)
+ (type bignum-element-type last))
+ (- digit-size (integer-length last) 1)))
+ ;;; Stores two bignums into the truncation bignum buffers,
+ ;;; shifting them on the way in. This assumes x and y are
+ ;;; positive and at least two in length, and it assumes
+ ;;; truncate-x and truncate-y are one digit longer than x and
+ ;;; y.
+ (shift-and-store-truncate-buffers (x len-x y len-y shift)
+ (declare (type bignum-index len-x len-y)
+ (type (integer 0 (#.digit-size)) shift))
+ (cond ((zerop shift)
+ (bignum-replace truncate-x x :end1 len-x)
+ (bignum-replace truncate-y y :end1 len-y))
+ (t
+ (bignum-ashift-left-unaligned x 0 shift (1+ len-x)
+ truncate-x)
+ (bignum-ashift-left-unaligned y 0 shift (1+ len-y)
+ truncate-y))))) ;; LABELS
+ ;;; Divide X by Y returning the quotient and remainder. In the
+ ;;; general case, we shift Y to set up for the algorithm, and we
+ ;;; use two buffers to save consing intermediate values. X gets
+ ;;; destructively modified to become the remainder, and we have
+ ;;; to shift it to account for the initial Y shift. After we
+ ;;; multiple bind q and r, we first fix up the signs and then
+ ;;; return the normalized results.
+ (let* ((x-plusp (%bignum-0-or-plusp x (%bignum-length x)))
+ (y-plusp (%bignum-0-or-plusp y (%bignum-length y)))
+ (x (if x-plusp x (negate-bignum x nil)))
+ (y (if y-plusp y (negate-bignum y nil)))
+ (len-x (%bignum-length x))
+ (len-y (%bignum-length y)))
+ (multiple-value-bind (q r)
+ (cond ((< len-y 2)
+ (bignum-truncate-single-digit x len-x y))
+ ((plusp (bignum-compare y x))
+ (let ((res (%allocate-bignum len-x)))
+ (dotimes (i len-x)
+ (setf (%bignum-ref res i) (%bignum-ref x i)))
+ (values 0 res)))
+ (t
+ (let ((len-x+1 (1+ len-x)))
+ (setf truncate-x (%allocate-bignum len-x+1))
+ (setf truncate-y (%allocate-bignum (1+ len-y)))
+ (let ((y-shift (shift-y-for-truncate y)))
+ (shift-and-store-truncate-buffers x len-x y
+ len-y y-shift)
+ (values (return-quotient-leaving-remainder len-x+1
+ len-y)
+ ;; Now that RETURN-QUOTIENT-LEAVING-REMAINDER
+ ;; has executed, we just tidy up the remainder
+ ;; (in TRUNCATE-X) and return it.
+ (cond
+ ((zerop y-shift)
+ (let ((res (%allocate-bignum len-y)))
+ (declare (type bignum-type res))
+ (bignum-replace res truncate-x :end2 len-y)
+ (%normalize-bignum res len-y)))
+ (t
+ (shift-right-unaligned
+ truncate-x 0 y-shift len-y
+ ((= j res-len-1)
+ (setf (%bignum-ref res j)
+ (%ashr (%bignum-ref truncate-x i)
+ y-shift))
+ (%normalize-bignum res res-len))
+ res))))))))
+ (let ((quotient (cond ((eq x-plusp y-plusp) q)
+ ((typep q 'fixnum) (the fixnum (- q)))
+ (t (negate-bignum-in-place q))))
+ (rem (cond (x-plusp r)
+ ((typep r 'fixnum) (the fixnum (- r)))
+ (t (negate-bignum-in-place r)))))
+ (values (if (typep quotient 'fixnum)
+ quotient
+ (%normalize-bignum quotient (%bignum-length quotient)))
+ (if (typep rem 'fixnum)
+ rem
+ (%normalize-bignum rem (%bignum-length rem))))))))))
-;;; This returns a guess for the next division step. Y1 is the highest y
-;;; digit, and y2 is the second to highest y digit. The x... variables are
-;;; the three highest x digits for the next division step.
-;;;
-;;; From Knuth, our guess is either all ones or x-i and x-i-1 divided by y1,
-;;; depending on whether x-i and y1 are the same. We test this guess by
-;;; determining whether guess*y2 is greater than the three high digits of x
-;;; minus guess*y1 shifted left one digit:
-;;; ------------------------------
-;;; | x-i | x-i-1 | x-i-2 |
-;;; ------------------------------
-;;; ------------------------------
-;;; - | g*y1 high | g*y1 low | 0 |
-;;; ------------------------------
-;;; ... < guess*y2 ???
-;;; If guess*y2 is greater, then we decrement our guess by one and try again.
-;;; This returns a guess that is either correct or one too large.
-(defun bignum-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
- (declare (type bignum-element-type y1 y2 x-i x-i-1 x-i-2))
- (let ((guess (if (%digit-compare x-i y1)
- all-ones-digit
- (%floor x-i x-i-1 y1))))
- (declare (type bignum-element-type guess))
- (loop
- (multiple-value-bind (high-guess*y1 low-guess*y1) (%multiply guess y1)
- (declare (type bignum-element-type low-guess*y1 high-guess*y1))
- (multiple-value-bind (high-guess*y2 low-guess*y2)
- (%multiply guess y2)
- (declare (type bignum-element-type high-guess*y2 low-guess*y2))
- (multiple-value-bind (middle-digit borrow)
- (%subtract-with-borrow x-i-1 low-guess*y1 1)
- (declare (type bignum-element-type middle-digit)
- (fixnum borrow))
- ;; Supplying borrow of 1 means there was no borrow, and we know
- ;; x-i-2 minus 0 requires no borrow.
- (let ((high-digit (%subtract-with-borrow x-i high-guess*y1 borrow)))
- (declare (type bignum-element-type high-digit))
- (if (and (%digit-compare high-digit 0)
- (or (%digit-greater high-guess*y2 middle-digit)
- (and (%digit-compare middle-digit high-guess*y2)
- (%digit-greater low-guess*y2 x-i-2))))
- (setf guess (%subtract-with-borrow guess 1 1))
- (return guess)))))))))
-
-;;; This returns the amount to shift y to place a one in the second highest
-;;; bit. Y must be positive. If the last digit of y is zero, then y has a
-;;; one in the previous digit's sign bit, so we know it will take one less
-;;; than digit-size to get a one where we want. Otherwise, we count how many
-;;; right shifts it takes to get zero; subtracting this value from digit-size
-;;; tells us how many high zeros there are which is one more than the shift
-;;; amount sought.
-;;;
-;;; Note: This is exactly the same as one less than the integer-length of the
-;;; last digit subtracted from the digit-size.
-;;;
-;;; We shift y to make it sufficiently large that doing the 64-bit by 32-bit
-;;; %FLOOR calls ensures the quotient and remainder fit in 32-bits.
-(defun shift-y-for-truncate (y)
- (let* ((len (%bignum-length y))
- (last (%bignum-ref y (1- len))))
- (declare (type bignum-index len)
- (type bignum-element-type last))
- (- digit-size (integer-length last) 1)))
-
-;;; Stores two bignums into the truncation bignum buffers, shifting them on the
-;;; way in. This assumes x and y are positive and at least two in length, and
-;;; it assumes *truncate-x* and *truncate-y* are one digit longer than x and y.
-(defun shift-and-store-truncate-buffers (x len-x y len-y shift)
- (declare (type bignum-index len-x len-y)
- (type (integer 0 (#.digit-size)) shift))
- (cond ((zerop shift)
- (bignum-replace *truncate-x* x :end1 len-x)
- (bignum-replace *truncate-y* y :end1 len-y))
- (t
- (bignum-ashift-left-unaligned x 0 shift (1+ len-x) *truncate-x*)
- (bignum-ashift-left-unaligned y 0 shift (1+ len-y) *truncate-y*))))
\f
;;;; %FLOOR primitive for BIGNUM-TRUNCATE
-;;; When a machine leaves out a 64-bit by 32-bit divide instruction (that is,
-;;; two bignum-digits divided by one), we have to roll our own (the hard way).
-;;; Basically, we treat the operation as four 16-bit digits divided by two
-;;; 16-bit digits. This means we have duplicated most of the code above to do
-;;; this nearly general 16-bit digit bignum divide, but we've unrolled loops
+;;; When a machine leaves out a 2*digit-size by digit-size divide
+;;; instruction (that is, two bignum-digits divided by one), we have to
+;;; roll our own (the hard way). Basically, we treat the operation as
+;;; four digit-size/2 digits divided by two digit-size/2 digits. This
+;;; means we have duplicated most of the code above to do this nearly
+;;; general digit-size/2 digit bignum divide, but we've unrolled loops
;;; and made use of other properties of this specific divide situation.
;;;; %FLOOR for machines with a 32x32 divider.
32x16-divide 32x16-multiply 32x16-multiply-split))
#!+32x16-divide
-(defconstant 32x16-base-1 #xFFFF)
+(defconstant 32x16-base-1 (1- (ash 1 (/ sb!vm:n-word-bits 2))))
-;;; This is similar to %SUBTRACT-WITH-BORROW. It returns a 16-bit difference
-;;; and a borrow. Returning a 1 for the borrow means there was no borrow, and
-;;; 0 means there was one.
+#!+32x16-divide
+(deftype bignum-half-element-type () `(unsigned-byte ,(/ sb!vm:n-word-bits 2)))
+#!+32x16-divide
+(defconstant half-digit-size (/ digit-size 2))
+
+;;; This is similar to %SUBTRACT-WITH-BORROW. It returns a
+;;; half-digit-size difference and a borrow. Returning a 1 for the
+;;; borrow means there was no borrow, and 0 means there was one.
#!+32x16-divide
(defun 32x16-subtract-with-borrow (a b borrow)
- (declare (type (unsigned-byte 16) a b)
+ (declare (type bignum-half-element-type a b)
(type (integer 0 1) borrow))
(let ((diff (+ (- a b) borrow 32x16-base-1)))
- (declare (type (unsigned-byte 17) diff))
- (values (logand diff #xFFFF)
- (ash diff -16))))
+ (declare (type (unsigned-byte #.(1+ half-digit-size)) diff))
+ (values (logand diff (1- (ash 1 half-digit-size)))
+ (ash diff (- half-digit-size)))))
-;;; This adds a and b, 16-bit quantities, with the carry k. It returns a
-;;; 16-bit sum and a second value, 0 or 1, indicating whether there was a
-;;; carry.
+;;; This adds a and b, half-digit-size quantities, with the carry k. It
+;;; returns a half-digit-size sum and a second value, 0 or 1, indicating
+;;; whether there was a carry.
#!+32x16-divide
(defun 32x16-add-with-carry (a b k)
- (declare (type (unsigned-byte 16) a b)
+ (declare (type bignum-half-element-type a b)
(type (integer 0 1) k))
(let ((res (the fixnum (+ a b k))))
- (declare (type (unsigned-byte 17) res))
- (if (zerop (the fixnum (logand #x10000 res)))
+ (declare (type (unsigned-byte #.(1+ half-digit-size)) res))
+ (if (zerop (the fixnum (logand (ash 1 half-digit-size) res)))
(values res 0)
- (values (the (unsigned-byte 16) (logand #xFFFF res))
+ (values (the bignum-half-element-type (logand (1- (ash 1 half-digit-size)) res))
1))))
-;;; This is probably a 32-bit by 32-bit divide instruction.
+;;; This is probably a digit-size by digit-size divide instruction.
#!+32x16-divide
(defun 32x16-divide (a b c)
- (declare (type (unsigned-byte 16) a b c))
+ (declare (type bignum-half-element-type a b c))
(floor (the bignum-element-type
(logior (the bignum-element-type (ash a 16))
b))
;;; register.
#!+32x16-divide
(defun 32x16-multiply (a b)
- (declare (type (unsigned-byte 16) a b))
+ (declare (type bignum-half-element-type a b))
(the bignum-element-type (* a b)))
-;;; This multiplies a and b, 16-bit quantities, and returns the result as two
-;;; 16-bit quantities, high and low.
+;;; This multiplies a and b, half-digit-size quantities, and returns the
+;;; result as two half-digit-size quantities, high and low.
#!+32x16-divide
(defun 32x16-multiply-split (a b)
(let ((res (32x16-multiply a b)))
(declare (the bignum-element-type res))
- (values (the (unsigned-byte 16) (logand #xFFFF (ash res -16)))
- (the (unsigned-byte 16) (logand #xFFFF res)))))
+ (values (the bignum-half-element-type (logand (1- (ash 1 half-digit-size)) (ash res (- half-digit-size))))
+ (the bignum-half-element-type (logand (1- (ash 1 half-digit-size)) res)))))
;;; The %FLOOR below uses this buffer the same way BIGNUM-TRUNCATE uses
-;;; *truncate-x*. There's no y buffer since we pass around the two 16-bit
-;;; digits and use them slightly differently than the general truncation
-;;; algorithm above.
+;;; *truncate-x*. There's no y buffer since we pass around the two
+;;; half-digit-size digits and use them slightly differently than the
+;;; general truncation algorithm above.
#!+32x16-divide
-(defvar *32x16-truncate-x* (make-array 4 :element-type '(unsigned-byte 16)
+(defvar *32x16-truncate-x* (make-array 4 :element-type 'bignum-half-element-type
:initial-element 0))
;;; This does the same thing as the %FLOOR above, but it does it at Lisp level
;;; when there is no 64x32-bit divide instruction on the machine.
;;;
-;;; It implements the higher level tactics of BIGNUM-TRUNCATE, but it makes use
-;;; of special situation provided, four 16-bit digits divided by two 16-bit
-;;; digits.
+;;; It implements the higher level tactics of BIGNUM-TRUNCATE, but it
+;;; makes use of special situation provided, four half-digit-size digits
+;;; divided by two half-digit-size digits.
#!+32x16-divide
(defun %floor (a b c)
(declare (type bignum-element-type a b c))
;; Setup *32x16-truncate-x* buffer from a and b.
(setf (aref *32x16-truncate-x* 0)
- (the (unsigned-byte 16) (logand #xFFFF b)))
+ (the bignum-half-element-type (logand (1- (ash 1 half-digit-size)) b)))
(setf (aref *32x16-truncate-x* 1)
- (the (unsigned-byte 16)
- (logand #xFFFF
- (the (unsigned-byte 16) (ash b -16)))))
+ (the bignum-half-element-type
+ (logand (1- (ash 1 half-digit-size))
+ (the bignum-half-element-type (ash b (- half-digit-size))))))
(setf (aref *32x16-truncate-x* 2)
- (the (unsigned-byte 16) (logand #xFFFF a)))
+ (the bignum-half-element-type (logand (1- (ash 1 half-digit-size)) a)))
(setf (aref *32x16-truncate-x* 3)
- (the (unsigned-byte 16)
- (logand #xFFFF
- (the (unsigned-byte 16) (ash a -16)))))
+ (the bignum-half-element-type
+ (logand (1- (ash 1 half-digit-size))
+ (the bignum-half-element-type (ash a (- half-digit-size))))))
;; From DO-TRUNCATE, but unroll the loop.
- (let* ((y1 (logand #xFFFF (ash c -16)))
- (y2 (logand #xFFFF c))
+ (let* ((y1 (logand (1- (ash 1 half-digit-size)) (ash c (- half-digit-size))))
+ (y2 (logand (1- (ash 1 half-digit-size)) c))
(q (the bignum-element-type
(ash (32x16-try-bignum-truncate-guess
(32x16-truncate-guess y1 y2
y1 y2 1)
16))))
(declare (type bignum-element-type q)
- (type (unsigned-byte 16) y1 y2))
+ (type bignum-half-element-type y1 y2))
(values (the bignum-element-type
(logior q
- (the (unsigned-byte 16)
+ (the bignum-half-element-type
(32x16-try-bignum-truncate-guess
(32x16-truncate-guess
y1 y2
(the bignum-element-type
(logior (the bignum-element-type
(ash (aref *32x16-truncate-x* 1) 16))
- (the (unsigned-byte 16)
+ (the bignum-half-element-type
(aref *32x16-truncate-x* 0)))))))
;;; This is similar to TRY-BIGNUM-TRUNCATE-GUESS, but this unrolls the two
#!+32x16-divide
(defun 32x16-try-bignum-truncate-guess (guess y-high y-low low-x-digit)
(declare (type bignum-index low-x-digit)
- (type (unsigned-byte 16) guess y-high y-low))
+ (type bignum-half-element-type guess y-high y-low))
(let ((high-x-digit (+ 2 low-x-digit)))
;; Multiply guess and divisor, subtracting from dividend simultaneously.
(multiple-value-bind (guess*y-hold carry borrow)
(32x16-try-guess-one-result-digit guess y-low 0 0 1 low-x-digit)
- (declare (type (unsigned-byte 16) guess*y-hold)
+ (declare (type bignum-half-element-type guess*y-hold)
(fixnum carry borrow))
(multiple-value-bind (guess*y-hold carry borrow)
(32x16-try-guess-one-result-digit guess y-high guess*y-hold
carry borrow (1+ low-x-digit))
- (declare (type (unsigned-byte 16) guess*y-hold)
+ (declare (type bignum-half-element-type guess*y-hold)
(fixnum borrow)
(ignore carry))
(setf (aref *32x16-truncate-x* high-x-digit)
(32x16-subtract-with-borrow (aref *32x16-truncate-x* high-x-digit)
guess*y-hold borrow))))
;; See whether guess is off by one, adding one Y back in if necessary.
- (cond ((zerop (logand #x8000 (aref *32x16-truncate-x* high-x-digit)))
+ (cond ((zerop (logand (ash 1 (1- half-digit-size))
+ (aref *32x16-truncate-x* high-x-digit)))
;; The subtraction result is zero or positive.
guess)
(t
(32x16-add-with-carry y-low
(aref *32x16-truncate-x* low-x-digit)
0)
- (declare (type (unsigned-byte 16) v))
+ (declare (type bignum-half-element-type v))
(setf (aref *32x16-truncate-x* low-x-digit) v)
(multiple-value-bind (v carry)
(32x16-add-with-carry y-high
(setf (aref *32x16-truncate-x* high-x-digit)
(32x16-add-with-carry (aref *32x16-truncate-x* high-x-digit)
carry 0))))
- (if (zerop (logand #x8000 guess))
+ (if (zerop (logand (ash 1 (1- half-digit-size)) guess))
(1- guess)
(1+ guess))))))
carry borrow x-index)
(multiple-value-bind (high-digit low-digit)
(32x16-multiply-split guess y-digit)
- (declare (type (unsigned-byte 16) high-digit low-digit))
+ (declare (type bignum-half-element-type high-digit low-digit))
(multiple-value-bind (low-digit temp-carry)
(32x16-add-with-carry low-digit guess*y-hold carry)
- (declare (type (unsigned-byte 16) low-digit))
+ (declare (type bignum-half-element-type low-digit))
(multiple-value-bind (high-digit temp-carry)
(32x16-add-with-carry high-digit temp-carry 0)
- (declare (type (unsigned-byte 16) high-digit))
+ (declare (type bignum-half-element-type high-digit))
(multiple-value-bind (x temp-borrow)
(32x16-subtract-with-borrow (aref *32x16-truncate-x* x-index)
low-digit borrow)
- (declare (type (unsigned-byte 16) x))
+ (declare (type bignum-half-element-type x))
(setf (aref *32x16-truncate-x* x-index) x)
(values high-digit temp-carry temp-borrow))))))
-;;; This is similar to BIGNUM-TRUNCATE-GUESS, but instead of computing the
-;;; guess exactly as described in the its comments (digit by digit), this
-;;; massages the 16-bit quantities into 32-bit quantities and performs the
+;;; This is similar to BIGNUM-TRUNCATE-GUESS, but instead of computing
+;;; the guess exactly as described in the its comments (digit by digit),
+;;; this massages the digit-size/2 quantities into digit-size quantities
+;;; and performs the
#!+32x16-divide
(defun 32x16-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
- (declare (type (unsigned-byte 16) y1 y2 x-i x-i-1 x-i-2))
+ (declare (type bignum-half-element-type y1 y2 x-i x-i-1 x-i-2))
(let ((guess (if (= x-i y1)
- #xFFFF
+ (1- (ash 1 half-digit-size))
(32x16-divide x-i x-i-1 y1))))
- (declare (type (unsigned-byte 16) guess))
+ (declare (type bignum-half-element-type guess))
(loop
(let* ((guess*y1 (the bignum-element-type
- (ash (logand #xFFFF
+ (ash (logand (1- (ash 1 half-digit-size))
(the bignum-element-type
(32x16-multiply guess y1)))
16)))
\f
;;;; general utilities
-;;; Allocate a single word bignum that holds fixnum. This is useful when
-;;; we are trying to mix fixnum and bignum operands.
-#!-sb-fluid (declaim (inline make-small-bignum))
-(defun make-small-bignum (fixnum)
- (let ((res (%allocate-bignum 1)))
- (setf (%bignum-ref res 0) (%fixnum-to-digit fixnum))
- res))
-
;;; Internal in-place operations use this to fixup remaining digits in the
;;; incoming data, such as in-place shifting. This is basically the same as
;;; the first form in %NORMALIZE-BIGNUM, but we return the length of the buffer
(%bignum-set-length result newlen))
(if (= newlen 1)
(let ((digit (%bignum-ref result 0)))
- (if (= (%ashr digit 29) (%ashr digit (1- digit-size)))
+ (if (= (%ashr digit sb!vm:n-positive-fixnum-bits)
+ (%ashr digit (1- digit-size)))
(%fixnum-digit-with-correct-sign digit)
result))
result)))
projects / sbcl.git / blobdiff