+;;; 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))
+