1 ;;;; code to implement bignum support
3 ;;;; This software is part of the SBCL system. See the README file for
6 ;;;; This software is derived from the CMU CL system, which was
7 ;;;; written at Carnegie Mellon University and released into the
8 ;;;; public domain. The software is in the public domain and is
9 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
10 ;;;; files for more information.
12 (in-package "SB!BIGNUM")
16 ;;; comments from CMU CL:
17 ;;; These symbols define the interface to the number code:
18 ;;; add-bignums multiply-bignums negate-bignum subtract-bignum
19 ;;; multiply-bignum-and-fixnum multiply-fixnums
20 ;;; bignum-ashift-right bignum-ashift-left bignum-gcd
21 ;;; bignum-to-float bignum-integer-length
22 ;;; bignum-logical-and bignum-logical-ior bignum-logical-xor
23 ;;; bignum-logical-not bignum-load-byte bignum-deposit-byte
24 ;;; bignum-truncate bignum-plus-p bignum-compare make-small-bignum
25 ;;; bignum-logbitp bignum-logcount
26 ;;; These symbols define the interface to the compiler:
27 ;;; bignum-type bignum-element-type bignum-index %allocate-bignum
28 ;;; %bignum-length %bignum-set-length %bignum-ref %bignum-set
29 ;;; %digit-0-or-plusp %add-with-carry %subtract-with-borrow
30 ;;; %multiply-and-add %multiply %lognot %logand %logior %logxor
31 ;;; %fixnum-to-digit %floor %fixnum-digit-with-correct-sign %ashl
32 ;;; %ashr %digit-logical-shift-right))
34 ;;; The following interfaces will either be assembler routines or code
35 ;;; sequences expanded into the code as basic bignum operations:
41 ;;; %BIGNUM-SET-LENGTH
42 ;;; %FIXNUM-DIGIT-WITH-CORRECT-SIGN
46 ;;; %BIGNUM-0-OR-PLUSP
47 ;;; %DIGIT-LOGICAL-SHIFT-RIGHT
48 ;;; General (May not exist when done due to sole use in %-routines.)
53 ;;; %SUBTRACT-WITH-BORROW
58 ;;; Shifting (in place)
59 ;;; %NORMALIZE-BIGNUM-BUFFER
60 ;;; GCD/Relational operators:
63 ;;; Relational operators:
72 ;;; Note: The floating routines know about the float representation.
75 ;;; There might be a problem with various LET's and parameters that take a
76 ;;; digit value. We need to write these so those things stay in machine
77 ;;; registers and number stack slots. I bind locals to these values, and I
78 ;;; use function on them -- ZEROP, ASH, etc.
81 ;;; In shifting and byte operations, I use masks and logical operations that
82 ;;; could result in intermediate bignums. This is hidden by the current system,
83 ;;; but I may need to write these in a way that keeps these masks and logical
84 ;;; operations from diving into the Lisp level bignum code.
88 ;;; logior, logxor, logand
89 ;;; depending on relationals, < (twice) and <= (twice)
90 ;;; or write compare thing (twice).
91 ;;; LDB on fixnum with bignum result.
92 ;;; DPB on fixnum with bignum result.
93 ;;; TRUNCATE returns zero or one as one value and fixnum or minus fixnum
94 ;;; for the other value when given (truncate fixnum bignum).
95 ;;; Returns (truncate bignum fixnum) otherwise.
97 ;;; subtraction (twice)
100 ;;; Write MASK-FIELD and DEPOSIT-FIELD in terms of logical operations.
102 ;;; IF (/ x y) with bignums:
103 ;;; do the truncate, and if rem is 0, return quotient.
106 ;;; "truncate" each by gcd, ignoring remainder 0.
107 ;;; form ratio of each result, bottom is positive.
109 ;;;; What's a bignum?
111 (defconstant digit-size sb!vm:n-word-bits)
113 (defconstant maximum-bignum-length (1- (ash 1 (- sb!vm:n-word-bits
114 sb!vm:n-widetag-bits))))
116 (defconstant all-ones-digit (1- (ash 1 sb!vm:n-word-bits)))
118 ;;;; internal inline routines
120 ;;; %ALLOCATE-BIGNUM must zero all elements.
121 (defun %allocate-bignum (length)
122 (declare (type bignum-index length))
123 (%allocate-bignum length))
125 ;;; Extract the length of the bignum.
126 (defun %bignum-length (bignum)
127 (declare (type bignum-type bignum))
128 (%bignum-length bignum))
130 ;;; %BIGNUM-REF needs to access bignums as obviously as possible, and it needs
131 ;;; to be able to return the digit somewhere no one looks for real objects.
132 (defun %bignum-ref (bignum i)
133 (declare (type bignum-type bignum)
134 (type bignum-index i))
135 (%bignum-ref bignum i))
136 (defun %bignum-set (bignum i value)
137 (declare (type bignum-type bignum)
138 (type bignum-index i)
139 (type bignum-element-type value))
140 (%bignum-set bignum i value))
142 ;;; Return T if digit is positive, or NIL if negative.
143 (defun %digit-0-or-plusp (digit)
144 (declare (type bignum-element-type digit))
145 (not (logbitp (1- digit-size) digit)))
147 #!-sb-fluid (declaim (inline %bignum-0-or-plusp))
148 (defun %bignum-0-or-plusp (bignum len)
149 (declare (type bignum-type bignum)
150 (type bignum-index len))
151 (%digit-0-or-plusp (%bignum-ref bignum (1- len))))
153 ;;; This should be in assembler, and should not cons intermediate
154 ;;; results. It returns a bignum digit and a carry resulting from adding
155 ;;; together a, b, and an incoming carry.
156 (defun %add-with-carry (a b carry)
157 (declare (type bignum-element-type a b)
158 (type (mod 2) carry))
159 (%add-with-carry a b carry))
161 ;;; This should be in assembler, and should not cons intermediate
162 ;;; results. It returns a bignum digit and a borrow resulting from
163 ;;; subtracting b from a, and subtracting a possible incoming borrow.
165 ;;; We really do: a - b - 1 + borrow, where borrow is either 0 or 1.
166 (defun %subtract-with-borrow (a b borrow)
167 (declare (type bignum-element-type a b)
168 (type (mod 2) borrow))
169 (%subtract-with-borrow a b borrow))
171 ;;; Multiply two digit-size numbers, returning a 2*digit-size result
172 ;;; split into two digit-size quantities.
173 (defun %multiply (x y)
174 (declare (type bignum-element-type x y))
177 ;;; This multiplies x-digit and y-digit, producing high and low digits
178 ;;; manifesting the result. Then it adds the low digit, res-digit, and
179 ;;; carry-in-digit. Any carries (note, you still have to add two digits
180 ;;; at a time possibly producing two carries) from adding these three
181 ;;; digits get added to the high digit from the multiply, producing the
182 ;;; next carry digit. Res-digit is optional since two uses of this
183 ;;; primitive multiplies a single digit bignum by a multiple digit
184 ;;; bignum, and in this situation there is no need for a result buffer
185 ;;; accumulating partial results which is where the res-digit comes
187 (defun %multiply-and-add (x-digit y-digit carry-in-digit
188 &optional (res-digit 0))
189 (declare (type bignum-element-type x-digit y-digit res-digit carry-in-digit))
190 (%multiply-and-add x-digit y-digit carry-in-digit res-digit))
192 (defun %lognot (digit)
193 (declare (type bignum-element-type digit))
196 ;;; Each of these does the digit-size unsigned op.
197 (declaim (inline %logand %logior %logxor))
199 (declare (type bignum-element-type a b))
202 (declare (type bignum-element-type a b))
205 (declare (type bignum-element-type a b))
208 ;;; This takes a fixnum and sets it up as an unsigned digit-size
210 (defun %fixnum-to-digit (x)
212 (logand x (1- (ash 1 digit-size))))
215 ;;; This takes three digits and returns the FLOOR'ed result of
216 ;;; dividing the first two as a 2*digit-size integer by the third.
218 ;;; Do weird LET and SETQ stuff to bamboozle the compiler into allowing
219 ;;; the %FLOOR transform to expand into pseudo-assembler for which the
220 ;;; compiler can later correctly allocate registers.
221 (defun %floor (a b c)
222 (let ((a a) (b b) (c c))
223 (declare (type bignum-element-type a b c))
227 ;;; Convert the digit to a regular integer assuming that the digit is signed.
228 (defun %fixnum-digit-with-correct-sign (digit)
229 (declare (type bignum-element-type digit))
230 (if (logbitp (1- digit-size) digit)
231 (logior digit (ash -1 digit-size))
234 ;;; Do an arithmetic shift right of data even though bignum-element-type is
236 (defun %ashr (data count)
237 (declare (type bignum-element-type data)
238 (type (mod #.sb!vm:n-word-bits) count))
241 ;;; This takes a digit-size quantity and shifts it to the left,
242 ;;; returning a digit-size quantity.
243 (defun %ashl (data count)
244 (declare (type bignum-element-type data)
245 (type (mod #.sb!vm:n-word-bits) count))
248 ;;; Do an unsigned (logical) right shift of a digit by Count.
249 (defun %digit-logical-shift-right (data count)
250 (declare (type bignum-element-type data)
251 (type (mod #.sb!vm:n-word-bits) count))
252 (%digit-logical-shift-right data count))
254 ;;; Change the length of bignum to be newlen. Newlen must be the same or
255 ;;; smaller than the old length, and any elements beyond newlen must be zeroed.
256 (defun %bignum-set-length (bignum newlen)
257 (declare (type bignum-type bignum)
258 (type bignum-index newlen))
259 (%bignum-set-length bignum newlen))
261 ;;; This returns 0 or "-1" depending on whether the bignum is positive. This
262 ;;; is suitable for infinite sign extension to complete additions,
263 ;;; subtractions, negations, etc. This cannot return a -1 represented as
264 ;;; a negative fixnum since it would then have to low zeros.
265 #!-sb-fluid (declaim (inline %sign-digit))
266 (defun %sign-digit (bignum len)
267 (declare (type bignum-type bignum)
268 (type bignum-index len))
269 (%ashr (%bignum-ref bignum (1- len)) (1- digit-size)))
271 ;;; These take two digit-size quantities and compare or contrast them
272 ;;; without wasting time with incorrect type checking.
273 (declaim (inline %digit-compare %digit-greater))
274 (defun %digit-compare (x y)
276 (defun %digit-greater (x y)
279 (declaim (optimize (speed 3) (safety 0)))
283 (defun add-bignums (a b)
284 (declare (type bignum-type a b))
285 (let ((len-a (%bignum-length a))
286 (len-b (%bignum-length b)))
287 (declare (type bignum-index len-a len-b))
288 (multiple-value-bind (a len-a b len-b)
290 (values a len-a b len-b)
291 (values b len-b a len-a))
292 (declare (type bignum-type a b)
293 (type bignum-index len-a len-b))
294 (let* ((len-res (1+ len-a))
295 (res (%allocate-bignum len-res))
297 (declare (type bignum-index len-res)
298 (type bignum-type res)
299 (type (mod 2) carry))
301 (declare (type bignum-index i))
302 (multiple-value-bind (v k)
303 (%add-with-carry (%bignum-ref a i) (%bignum-ref b i) carry)
304 (declare (type bignum-element-type v)
306 (setf (%bignum-ref res i) v)
309 (finish-add a res carry (%sign-digit b len-b) len-b len-a)
310 (setf (%bignum-ref res len-a)
311 (%add-with-carry (%sign-digit a len-a)
312 (%sign-digit b len-b)
314 (%normalize-bignum res len-res)))))
316 ;;; This takes the longer of two bignums and propagates the carry through its
317 ;;; remaining high order digits.
318 (defun finish-add (a res carry sign-digit-b start end)
319 (declare (type bignum-type a res)
321 (type bignum-element-type sign-digit-b)
322 (type bignum-index start end))
323 (do ((i start (1+ i)))
325 (setf (%bignum-ref res end)
326 (%add-with-carry (%sign-digit a end) sign-digit-b carry)))
327 (declare (type bignum-index i))
328 (multiple-value-bind (v k)
329 (%add-with-carry (%bignum-ref a i) sign-digit-b carry)
330 (setf (%bignum-ref res i) v)
336 (eval-when (:compile-toplevel :execute)
338 ;;; This subtracts b from a plugging result into res. Return-fun is the
339 ;;; function to call that fixes up the result returning any useful values, such
340 ;;; as the result. This macro may evaluate its arguments more than once.
341 (sb!xc:defmacro subtract-bignum-loop (a len-a b len-b res len-res return-fun)
342 (with-unique-names (borrow a-digit a-sign b-digit b-sign i v k)
344 (,a-sign (%sign-digit ,a ,len-a))
345 (,b-sign (%sign-digit ,b ,len-b)))
346 (declare (type bignum-element-type ,a-sign ,b-sign))
347 (dotimes (,i ,len-res)
348 (declare (type bignum-index ,i))
349 (let ((,a-digit (if (< ,i ,len-a) (%bignum-ref ,a ,i) ,a-sign))
350 (,b-digit (if (< ,i ,len-b) (%bignum-ref ,b ,i) ,b-sign)))
351 (declare (type bignum-element-type ,a-digit ,b-digit))
352 (multiple-value-bind (,v ,k)
353 (%subtract-with-borrow ,a-digit ,b-digit ,borrow)
354 (setf (%bignum-ref ,res ,i) ,v)
356 (,return-fun ,res ,len-res))))
360 (defun subtract-bignum (a b)
361 (declare (type bignum-type a b))
362 (let* ((len-a (%bignum-length a))
363 (len-b (%bignum-length b))
364 (len-res (1+ (max len-a len-b)))
365 (res (%allocate-bignum len-res)))
366 (declare (type bignum-index len-a len-b len-res)) ;Test len-res for bounds?
367 (subtract-bignum-loop a len-a b len-b res len-res %normalize-bignum)))
369 ;;; Operations requiring a subtraction without the overhead of intermediate
370 ;;; results, such as GCD, use this. It assumes Result is big enough for the
372 (defun subtract-bignum-buffers-with-len (a len-a b len-b result len-res)
373 (declare (type bignum-type a b result)
374 (type bignum-index len-a len-b len-res))
375 (subtract-bignum-loop a len-a b len-b result len-res
376 %normalize-bignum-buffer))
378 (defun subtract-bignum-buffers (a len-a b len-b result)
379 (declare (type bignum-type a b result)
380 (type bignum-index len-a len-b))
381 (subtract-bignum-loop a len-a b len-b result (max len-a len-b)
382 %normalize-bignum-buffer))
386 (defun multiply-bignums (a b)
387 (declare (type bignum-type a b))
388 (let* ((a-plusp (%bignum-0-or-plusp a (%bignum-length a)))
389 (b-plusp (%bignum-0-or-plusp b (%bignum-length b)))
390 (a (if a-plusp a (negate-bignum a)))
391 (b (if b-plusp b (negate-bignum b)))
392 (len-a (%bignum-length a))
393 (len-b (%bignum-length b))
394 (len-res (+ len-a len-b))
395 (res (%allocate-bignum len-res))
396 (negate-res (not (eq a-plusp b-plusp))))
397 (declare (type bignum-index len-a len-b len-res))
399 (declare (type bignum-index i))
400 (let ((carry-digit 0)
401 (x (%bignum-ref a i))
403 (declare (type bignum-index k)
404 (type bignum-element-type carry-digit x))
406 (multiple-value-bind (big-carry res-digit)
411 (declare (type bignum-element-type big-carry res-digit))
412 (setf (%bignum-ref res k) res-digit)
413 (setf carry-digit big-carry)
415 (setf (%bignum-ref res k) carry-digit)))
416 (when negate-res (negate-bignum-in-place res))
417 (%normalize-bignum res len-res)))
419 (defun multiply-bignum-and-fixnum (bignum fixnum)
420 (declare (type bignum-type bignum) (type fixnum fixnum))
421 (let* ((bignum-plus-p (%bignum-0-or-plusp bignum (%bignum-length bignum)))
422 (fixnum-plus-p (not (minusp fixnum)))
423 (bignum (if bignum-plus-p bignum (negate-bignum bignum)))
424 (bignum-len (%bignum-length bignum))
425 (fixnum (if fixnum-plus-p fixnum (- fixnum)))
426 (result (%allocate-bignum (1+ bignum-len)))
428 (declare (type bignum-type bignum result)
429 (type bignum-index bignum-len)
430 (type bignum-element-type fixnum carry-digit))
431 (dotimes (index bignum-len)
432 (declare (type bignum-index index))
433 (multiple-value-bind (next-digit low)
434 (%multiply-and-add (%bignum-ref bignum index) fixnum carry-digit)
435 (declare (type bignum-element-type next-digit low))
436 (setf carry-digit next-digit)
437 (setf (%bignum-ref result index) low)))
438 (setf (%bignum-ref result bignum-len) carry-digit)
439 (unless (eq bignum-plus-p fixnum-plus-p)
440 (negate-bignum-in-place result))
441 (%normalize-bignum result (1+ bignum-len))))
443 (defun multiply-fixnums (a b)
444 (declare (fixnum a b))
445 (let* ((a-minusp (minusp a))
446 (b-minusp (minusp b)))
447 (multiple-value-bind (high low)
448 (%multiply (if a-minusp (- a) a)
449 (if b-minusp (- b) b))
450 (declare (type bignum-element-type high low))
451 (if (and (zerop high)
452 (%digit-0-or-plusp low))
453 (let ((low (sb!ext:truly-the (unsigned-byte #.(1- sb!vm:n-word-bits))
454 (%fixnum-digit-with-correct-sign low))))
455 (if (eq a-minusp b-minusp)
458 (let ((res (%allocate-bignum 2)))
459 (%bignum-set res 0 low)
460 (%bignum-set res 1 high)
461 (unless (eq a-minusp b-minusp) (negate-bignum-in-place res))
462 (%normalize-bignum res 2))))))
464 ;;;; BIGNUM-REPLACE and WITH-BIGNUM-BUFFERS
466 (eval-when (:compile-toplevel :execute)
468 (sb!xc:defmacro bignum-replace (dest
476 (sb!int:once-only ((n-dest dest)
478 (with-unique-names (n-start1 n-end1 n-start2 n-end2 i1 i2)
479 (let ((end1 (or end1 `(%bignum-length ,n-dest)))
480 (end2 (or end2 `(%bignum-length ,n-src))))
482 `(let ((,n-start1 ,start1)
484 (do ((,i1 (1- ,end1) (1- ,i1))
485 (,i2 (1- ,end2) (1- ,i2)))
486 ((or (< ,i1 ,n-start1) (< ,i2 ,n-start2)))
487 (declare (fixnum ,i1 ,i2))
488 (%bignum-set ,n-dest ,i1 (%bignum-ref ,n-src ,i2))))
489 (if (eql start1 start2)
490 `(let ((,n-end1 (min ,end1 ,end2)))
491 (do ((,i1 ,start1 (1+ ,i1)))
493 (declare (type bignum-index ,i1))
494 (%bignum-set ,n-dest ,i1 (%bignum-ref ,n-src ,i1))))
495 `(let ((,n-end1 ,end1)
497 (do ((,i1 ,start1 (1+ ,i1))
498 (,i2 ,start2 (1+ ,i2)))
499 ((or (>= ,i1 ,n-end1) (>= ,i2 ,n-end2)))
500 (declare (type bignum-index ,i1 ,i2))
501 (%bignum-set ,n-dest ,i1 (%bignum-ref ,n-src ,i2))))))))))
503 (sb!xc:defmacro with-bignum-buffers (specs &body body)
505 "WITH-BIGNUM-BUFFERS ({(var size [init])}*) Form*"
506 (sb!int:collect ((binds)
509 (let ((name (first spec))
510 (size (second spec)))
511 (binds `(,name (%allocate-bignum ,size)))
512 (let ((init (third spec)))
514 (inits `(bignum-replace ,name ,init))))))
523 (eval-when (:compile-toplevel :load-toplevel :execute)
524 ;; The asserts in the GCD implementation are way too expensive to
525 ;; check in normal use, and are disabled here.
526 (sb!xc:defmacro gcd-assert (&rest args)
529 ;; We'll be doing a lot of modular arithmetic.
530 (sb!xc:defmacro modularly (form)
531 `(logand all-ones-digit ,form)))
533 ;;; I'm not sure why I need this FTYPE declaration. Compiled by the
534 ;;; target compiler, it can deduce the return type fine, but without
535 ;;; it, we pay a heavy price in BIGNUM-GCD when compiled by the
536 ;;; cross-compiler. -- CSR, 2004-07-19
537 (declaim (ftype (sfunction (bignum-type bignum-index bignum-type bignum-index)
538 (and unsigned-byte fixnum))
539 bignum-factors-of-two))
540 (defun bignum-factors-of-two (a len-a b len-b)
541 (declare (type bignum-index len-a len-b) (type bignum-type a b))
543 (end (min len-a len-b)))
544 ((= i end) (error "Unexpected zero bignums?"))
545 (declare (type bignum-index i end))
546 (let ((or-digits (%logior (%bignum-ref a i) (%bignum-ref b i))))
547 (unless (zerop or-digits)
548 (return (do ((j 0 (1+ j))
549 (or-digits or-digits (%ashr or-digits 1)))
550 ((oddp or-digits) (+ (* i digit-size) j))
551 (declare (type (mod #.sb!vm:n-word-bits) j))))))))
553 ;;; Multiply a bignum buffer with a fixnum or a digit, storing the
554 ;;; result in another bignum buffer, and without using any
555 ;;; temporaries. Inlined to avoid boxing smallnum if it's actually a
556 ;;; digit. Needed by GCD, should possibly OAOO with
557 ;;; MULTIPLY-BIGNUM-AND-FIXNUM.
558 (declaim (inline multiply-bignum-buffer-and-smallnum-to-buffer))
559 (defun multiply-bignum-buffer-and-smallnum-to-buffer (bignum bignum-len
561 (declare (type bignum-type bignum))
562 (let* ((bignum-plus-p (%bignum-0-or-plusp bignum bignum-len))
563 (smallnum-plus-p (not (minusp smallnum)))
564 (smallnum (if smallnum-plus-p smallnum (- smallnum)))
566 (declare (type bignum-type bignum res)
567 (type bignum-index bignum-len)
568 (type bignum-element-type smallnum carry-digit))
569 (unless bignum-plus-p
570 (negate-bignum-buffer-in-place bignum bignum-len))
571 (dotimes (index bignum-len)
572 (declare (type bignum-index index))
573 (multiple-value-bind (next-digit low)
574 (%multiply-and-add (%bignum-ref bignum index)
577 (declare (type bignum-element-type next-digit low))
578 (setf carry-digit next-digit)
579 (setf (%bignum-ref res index) low)))
580 (setf (%bignum-ref res bignum-len) carry-digit)
581 (unless bignum-plus-p
582 (negate-bignum-buffer-in-place bignum bignum-len))
583 (let ((res-len (%normalize-bignum-buffer res (1+ bignum-len))))
584 (unless (eq bignum-plus-p smallnum-plus-p)
585 (negate-bignum-buffer-in-place res res-len))
588 ;;; Given U and V, return U / V mod 2^32. Implements the algorithm in the
589 ;;; paper, but uses some clever bit-twiddling nicked from Nickle to do it.
590 (declaim (inline bmod))
592 (let ((ud (%bignum-ref u 0))
593 (vd (%bignum-ref v 0))
597 (declare (type (unsigned-byte #.sb!vm:n-word-bits) ud vd umask imask m))
598 (dotimes (i digit-size)
599 (setf umask (logior umask imask))
600 (when (logtest ud umask)
601 (setf ud (modularly (- ud vd)))
602 (setf m (modularly (logior m imask))))
603 (setf imask (modularly (ash imask 1)))
604 (setf vd (modularly (ash vd 1))))
607 (defun dmod (u u-len v v-len tmp1)
608 (loop while (> (bignum-buffer-integer-length u u-len)
609 (+ (bignum-buffer-integer-length v v-len)
612 (unless (zerop (%bignum-ref u 0))
613 (let* ((bmod (bmod u v))
614 (tmp1-len (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
617 (setf u-len (subtract-bignum-buffers u u-len
620 (bignum-abs-buffer u u-len)))
621 (gcd-assert (zerop (%bignum-ref u 0)))
622 (setf u-len (bignum-buffer-ashift-right u u-len digit-size)))
623 (let* ((d (+ 1 (- (bignum-buffer-integer-length u u-len)
624 (bignum-buffer-integer-length v v-len))))
626 (declare (type (unsigned-byte #.(integer-length #.sb!vm:n-word-bits)) d)
627 (type (unsigned-byte #.sb!vm:n-word-bits) n))
628 (gcd-assert (>= d 0))
629 (when (logtest (%bignum-ref u 0) n)
631 (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
635 (setf u-len (subtract-bignum-buffers u u-len
638 (bignum-abs-buffer u u-len)))
641 (defconstant lower-ones-digit (1- (ash 1 (truncate sb!vm:n-word-bits 2))))
643 ;;; Find D and N such that (LOGAND ALL-ONES-DIGIT (- (* D X) (* N Y))) is 0,
644 ;;; (< 0 N LOWER-ONES-DIGIT) and (< 0 (ABS D) LOWER-ONES-DIGIT).
645 (defun reduced-ratio-mod (x y)
646 (let* ((c (bmod x y))
649 (n2 (modularly (1+ (modularly (lognot n1)))))
651 (declare (type (unsigned-byte #.sb!vm:n-word-bits) n1 d1 n2 d2))
652 (loop while (> n2 (expt 2 (truncate digit-size 2))) do
653 (loop for i of-type (mod #.sb!vm:n-word-bits)
654 downfrom (- (integer-length n1) (integer-length n2))
656 (when (>= n1 (modularly (ash n2 i)))
657 (psetf n1 (modularly (- n1 (modularly (ash n2 i))))
658 d1 (modularly (- d1 (modularly (ash d2 i)))))))
663 (values n2 (if (>= d2 (expt 2 (1- digit-size)))
664 (lognot (logand most-positive-fixnum (lognot d2)))
665 (logand lower-ones-digit d2)))))
668 (defun copy-bignum (a &optional (len (%bignum-length a)))
669 (let ((b (%allocate-bignum len)))
671 (%bignum-set-length b len)
674 ;;; Allocate a single word bignum that holds fixnum. This is useful when
675 ;;; we are trying to mix fixnum and bignum operands.
676 #!-sb-fluid (declaim (inline make-small-bignum))
677 (defun make-small-bignum (fixnum)
678 (let ((res (%allocate-bignum 1)))
679 (setf (%bignum-ref res 0) (%fixnum-to-digit fixnum))
682 ;; When the larger number is less than this many bignum digits long, revert
684 (defparameter *accelerated-gcd-cutoff* 3)
686 ;;; Alternate between k-ary reduction with the help of
687 ;;; REDUCED-RATIO-MOD and digit modulus reduction via DMOD. Once the
688 ;;; arguments get small enough, drop through to BIGNUM-MOD-GCD (since
689 ;;; k-ary reduction can introduce spurious factors, which need to be
690 ;;; filtered out). Reference: Kenneth Weber, "The accelerated integer
691 ;;; GCD algorithm", ACM Transactions on Mathematical Software, volume
692 ;;; 21, number 1, March 1995, epp. 111-122.
693 (defun bignum-gcd (u0 v0)
694 (declare (type bignum-type u0 v0))
695 (let* ((u1 (if (%bignum-0-or-plusp u0 (%bignum-length u0))
697 (negate-bignum u0 nil)))
698 (v1 (if (%bignum-0-or-plusp v0 (%bignum-length v0))
700 (negate-bignum v0 nil))))
702 (return-from bignum-gcd u1))
705 (let ((n (mod v1 u1)))
706 (setf v1 (if (fixnump n)
707 (make-small-bignum n)
709 (if (and (= 1 (%bignum-length v1))
710 (zerop (%bignum-ref v1 0)))
711 (return-from bignum-gcd (%normalize-bignum u1
712 (%bignum-length u1))))
713 (let* ((buffer-len (+ 2 (%bignum-length u1)))
714 (u (%allocate-bignum buffer-len))
715 (u-len (%bignum-length u1))
716 (v (%allocate-bignum buffer-len))
717 (v-len (%bignum-length v1))
718 (tmp1 (%allocate-bignum buffer-len))
720 (tmp2 (%allocate-bignum buffer-len))
723 (bignum-factors-of-two u1 (%bignum-length u1)
724 v1 (%bignum-length v1))))
725 (declare (type (or null bignum-index)
726 buffer-len u-len v-len tmp1-len tmp2-len))
727 (bignum-replace u u1)
728 (bignum-replace v v1)
730 (make-gcd-bignum-odd u
731 (bignum-buffer-ashift-right u u-len
734 (make-gcd-bignum-odd v
735 (bignum-buffer-ashift-right v v-len
737 (loop until (or (< u-len *accelerated-gcd-cutoff*)
741 (zerop (%bignum-ref v 0))))
743 (gcd-assert (= buffer-len (%bignum-length u)
745 (%bignum-length tmp1)
746 (%bignum-length tmp2)))
747 (if (> (bignum-buffer-integer-length u u-len)
748 (+ #.(truncate sb!vm:n-word-bits 4)
749 (bignum-buffer-integer-length v v-len)))
750 (setf u-len (dmod u u-len
753 (multiple-value-bind (n d) (reduced-ratio-mod u v)
755 (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
758 (multiply-bignum-buffer-and-smallnum-to-buffer u u-len
760 (gcd-assert (= (copy-bignum tmp2 tmp2-len)
761 (* (copy-bignum u u-len) d)))
762 (gcd-assert (= (copy-bignum tmp1 tmp1-len)
763 (* (copy-bignum v v-len) n)))
765 (subtract-bignum-buffers-with-len tmp1 tmp1-len
770 (gcd-assert (or (zerop (- (copy-bignum tmp1 tmp1-len)
771 (copy-bignum tmp2 tmp2-len)))
772 (= (copy-bignum u u-len)
773 (- (copy-bignum tmp1 tmp1-len)
774 (copy-bignum tmp2 tmp2-len)))))
775 (bignum-abs-buffer u u-len)
776 (gcd-assert (zerop (modularly u)))))
777 (setf u-len (make-gcd-bignum-odd u u-len))
779 (rotatef u-len v-len))
780 (setf u (copy-bignum u u-len))
781 (let ((n (bignum-mod-gcd v1 u)))
782 (ash (bignum-mod-gcd u1 (if (fixnump n)
783 (make-small-bignum n)
787 (defun bignum-mod-gcd (a b)
788 (declare (type bignum-type a b))
791 ;; While the length difference of A and B is sufficiently large,
792 ;; reduce using MOD (slowish, but it should equalize the sizes of
793 ;; A and B pretty quickly). After that, use the binary GCD
794 ;; algorithm to handle the rest.
795 (loop until (and (= (%bignum-length b) 1) (zerop (%bignum-ref b 0))) do
796 (when (<= (%bignum-length a) (1+ (%bignum-length b)))
797 (return-from bignum-mod-gcd (bignum-binary-gcd a b)))
798 (let ((rem (mod a b)))
800 (setf a (make-small-bignum rem))
803 (if (= (%bignum-length a) 1)
804 (%normalize-bignum a 1)
807 (defun bignum-binary-gcd (a b)
808 (declare (type bignum-type a b))
809 (let* ((len-a (%bignum-length a))
810 (len-b (%bignum-length b)))
811 (declare (type bignum-index len-a len-b))
812 (with-bignum-buffers ((a-buffer len-a a)
814 (res-buffer (max len-a len-b)))
815 (let* ((factors-of-two
816 (bignum-factors-of-two a-buffer len-a
818 (len-a (make-gcd-bignum-odd
820 (bignum-buffer-ashift-right a-buffer len-a
822 (len-b (make-gcd-bignum-odd
824 (bignum-buffer-ashift-right b-buffer len-b
826 (declare (type bignum-index len-a len-b))
833 (multiple-value-bind (u v len-v r len-r)
834 (bignum-gcd-order-and-subtract x len-x y len-y z)
835 (declare (type bignum-index len-v len-r))
836 (when (and (= len-r 1) (zerop (%bignum-ref r 0)))
837 (if (zerop factors-of-two)
838 (let ((ret (%allocate-bignum len-v)))
840 (setf (%bignum-ref ret i) (%bignum-ref v i)))
841 (return (%normalize-bignum ret len-v)))
842 (return (bignum-ashift-left v factors-of-two len-v))))
843 (setf x v len-x len-v)
844 (setf y r len-y (make-gcd-bignum-odd r len-r))
847 (defun bignum-gcd-order-and-subtract (a len-a b len-b res)
848 (declare (type bignum-index len-a len-b) (type bignum-type a b))
849 (cond ((= len-a len-b)
850 (do ((i (1- len-a) (1- i)))
852 (setf (%bignum-ref res 0) 0)
853 (values a b len-b res 1))
854 (let ((a-digit (%bignum-ref a i))
855 (b-digit (%bignum-ref b i)))
856 (cond ((%digit-compare a-digit b-digit))
857 ((%digit-greater a-digit b-digit)
859 (values a b len-b res
860 (subtract-bignum-buffers a len-a b len-b
864 (values b a len-a res
865 (subtract-bignum-buffers b len-b
869 (values a b len-b res
870 (subtract-bignum-buffers a len-a b len-b res)))
872 (values b a len-a res
873 (subtract-bignum-buffers b len-b a len-a res)))))
875 (defun make-gcd-bignum-odd (a len-a)
876 (declare (type bignum-type a) (type bignum-index len-a))
877 (dotimes (index len-a)
878 (declare (type bignum-index index))
879 (do ((digit (%bignum-ref a index) (%ashr digit 1))
880 (increment 0 (1+ increment)))
882 (declare (type (mod #.sb!vm:n-word-bits) increment))
884 (return-from make-gcd-bignum-odd
885 (bignum-buffer-ashift-right a len-a
886 (+ (* index digit-size)
892 (eval-when (:compile-toplevel :execute)
894 ;;; This negates bignum-len digits of bignum, storing the resulting digits into
895 ;;; result (possibly EQ to bignum) and returning whatever end-carry there is.
896 (sb!xc:defmacro bignum-negate-loop
897 (bignum bignum-len &optional (result nil resultp))
898 (with-unique-names (carry end value last)
899 `(let* (,@(if (not resultp) `(,last))
901 (multiple-value-bind (,value ,carry)
902 (%add-with-carry (%lognot (%bignum-ref ,bignum 0)) 1 0)
904 `(setf (%bignum-ref ,result 0) ,value)
905 `(setf ,last ,value))
909 (declare (type bit ,carry)
910 (type bignum-index i ,end))
912 (when (= i ,end) (return))
913 (multiple-value-bind (,value temp)
914 (%add-with-carry (%lognot (%bignum-ref ,bignum i)) 0 ,carry)
916 `(setf (%bignum-ref ,result i) ,value)
917 `(setf ,last ,value))
920 ,(if resultp carry `(values ,carry ,last)))))
924 ;;; Fully-normalize is an internal optional. It cause this to always return
925 ;;; a bignum, without any extraneous digits, and it never returns a fixnum.
926 (defun negate-bignum (x &optional (fully-normalize t))
927 (declare (type bignum-type x))
928 (let* ((len-x (%bignum-length x))
930 (res (%allocate-bignum len-res)))
931 (declare (type bignum-index len-x len-res)) ;Test len-res for range?
932 (let ((carry (bignum-negate-loop x len-x res)))
933 (setf (%bignum-ref res len-x)
934 (%add-with-carry (%lognot (%sign-digit x len-x)) 0 carry)))
936 (%normalize-bignum res len-res)
937 (%mostly-normalize-bignum res len-res))))
939 ;;; This assumes bignum is positive; that is, the result of negating it will
940 ;;; stay in the provided allocated bignum.
941 (defun negate-bignum-buffer-in-place (bignum bignum-len)
942 (bignum-negate-loop bignum bignum-len bignum)
945 (defun negate-bignum-in-place (bignum)
946 (declare (inline negate-bignum-buffer-in-place))
947 (negate-bignum-buffer-in-place bignum (%bignum-length bignum)))
949 (defun bignum-abs-buffer (bignum len)
950 (unless (%bignum-0-or-plusp bignum len)
951 (negate-bignum-buffer-in-place bignum len)))
955 (eval-when (:compile-toplevel :execute)
957 ;;; This macro is used by BIGNUM-ASHIFT-RIGHT, BIGNUM-BUFFER-ASHIFT-RIGHT, and
958 ;;; BIGNUM-LDB-BIGNUM-RES. They supply a termination form that references
959 ;;; locals established by this form. Source is the source bignum. Start-digit
960 ;;; is the first digit in source from which we pull bits. Start-pos is the
961 ;;; first bit we want. Res-len-form is the form that computes the length of
962 ;;; the resulting bignum. Termination is a DO termination form with a test and
963 ;;; body. When result is supplied, it is the variable to which this binds a
964 ;;; newly allocated bignum.
966 ;;; Given start-pos, 1-31 inclusively, of shift, we form the j'th resulting
967 ;;; digit from high bits of the i'th source digit and the start-pos number of
968 ;;; bits from the i+1'th source digit.
969 (sb!xc:defmacro shift-right-unaligned (source
975 `(let* ((high-bits-in-first-digit (- digit-size ,start-pos))
976 (res-len ,res-len-form)
977 (res-len-1 (1- res-len))
978 ,@(if result `((,result (%allocate-bignum res-len)))))
979 (declare (type bignum-index res-len res-len-1))
980 (do ((i ,start-digit (1+ i))
983 (declare (type bignum-index i j))
984 (setf (%bignum-ref ,(if result result source) j)
985 (%logior (%digit-logical-shift-right (%bignum-ref ,source i)
987 (%ashl (%bignum-ref ,source (1+ i))
988 high-bits-in-first-digit))))))
992 ;;; First compute the number of whole digits to shift, shifting them by
993 ;;; skipping them when we start to pick up bits, and the number of bits to
994 ;;; shift the remaining digits into place. If the number of digits is greater
995 ;;; than the length of the bignum, then the result is either 0 or -1. If we
996 ;;; shift on a digit boundary (that is, n-bits is zero), then we just copy
997 ;;; digits. The last branch handles the general case which uses a macro that a
998 ;;; couple other routines use. The fifth argument to the macro references
999 ;;; locals established by the macro.
1000 (defun bignum-ashift-right (bignum count)
1001 (declare (type bignum-type bignum)
1002 (type unsigned-byte count))
1003 (let ((bignum-len (%bignum-length bignum)))
1004 (declare (type bignum-index bignum-len))
1005 (cond ((fixnump count)
1006 (multiple-value-bind (digits n-bits) (truncate count digit-size)
1007 (declare (type bignum-index digits))
1009 ((>= digits bignum-len)
1010 (if (%bignum-0-or-plusp bignum bignum-len) 0 -1))
1012 (bignum-ashift-right-digits bignum digits))
1014 (shift-right-unaligned bignum digits n-bits (- bignum-len digits)
1016 (setf (%bignum-ref res j)
1017 (%ashr (%bignum-ref bignum i) n-bits))
1018 (%normalize-bignum res res-len))
1020 ((> count bignum-len)
1021 (if (%bignum-0-or-plusp bignum bignum-len) 0 -1))
1022 ;; Since a FIXNUM should be big enough to address anything in
1023 ;; memory, including arrays of bits, and since arrays of bits
1024 ;; take up about the same space as corresponding fixnums, there
1025 ;; should be no way that we fall through to this case: any shift
1026 ;; right by a bignum should give zero. But let's check anyway:
1027 (t (error "bignum overflow: can't shift right by ~S" count)))))
1029 (defun bignum-ashift-right-digits (bignum digits)
1030 (declare (type bignum-type bignum)
1031 (type bignum-index digits))
1032 (let* ((res-len (- (%bignum-length bignum) digits))
1033 (res (%allocate-bignum res-len)))
1034 (declare (type bignum-index res-len)
1035 (type bignum-type res))
1036 (bignum-replace res bignum :start2 digits)
1037 (%normalize-bignum res res-len)))
1039 ;;; GCD uses this for an in-place shifting operation. This is different enough
1040 ;;; from BIGNUM-ASHIFT-RIGHT that it isn't worth folding the bodies into a
1041 ;;; macro, but they share the basic algorithm. This routine foregoes a first
1042 ;;; test for digits being greater than or equal to bignum-len since that will
1043 ;;; never happen for its uses in GCD. We did fold the last branch into a macro
1044 ;;; since it was duplicated a few times, and the fifth argument to it
1045 ;;; references locals established by the macro.
1046 (defun bignum-buffer-ashift-right (bignum bignum-len x)
1047 (declare (type bignum-index bignum-len) (fixnum x))
1048 (multiple-value-bind (digits n-bits) (truncate x digit-size)
1049 (declare (type bignum-index digits))
1052 (let ((new-end (- bignum-len digits)))
1053 (bignum-replace bignum bignum :end1 new-end :start2 digits
1055 (%normalize-bignum-buffer bignum new-end)))
1057 (shift-right-unaligned bignum digits n-bits (- bignum-len digits)
1059 (setf (%bignum-ref bignum j)
1060 (%ashr (%bignum-ref bignum i) n-bits))
1061 (%normalize-bignum-buffer bignum res-len)))))))
1063 ;;; This handles shifting a bignum buffer to provide fresh bignum data for some
1064 ;;; internal routines. We know bignum is safe when called with bignum-len.
1065 ;;; First we compute the number of whole digits to shift, shifting them
1066 ;;; starting to store farther along the result bignum. If we shift on a digit
1067 ;;; boundary (that is, n-bits is zero), then we just copy digits. The last
1068 ;;; branch handles the general case.
1069 (defun bignum-ashift-left (bignum x &optional bignum-len)
1070 (declare (type bignum-type bignum)
1071 (type unsigned-byte x)
1072 (type (or null bignum-index) bignum-len))
1074 (multiple-value-bind (digits n-bits) (truncate x digit-size)
1075 (let* ((bignum-len (or bignum-len (%bignum-length bignum)))
1076 (res-len (+ digits bignum-len 1)))
1077 (when (> res-len maximum-bignum-length)
1078 (error "can't represent result of left shift"))
1080 (bignum-ashift-left-digits bignum bignum-len digits)
1081 (bignum-ashift-left-unaligned bignum digits n-bits res-len))))
1082 ;; Left shift by a number too big to be represented as a fixnum
1083 ;; would exceed our memory capacity, since a fixnum is big enough
1084 ;; to index any array, including a bit array.
1085 (error "can't represent result of left shift")))
1087 (defun bignum-ashift-left-digits (bignum bignum-len digits)
1088 (declare (type bignum-index bignum-len digits))
1089 (let* ((res-len (+ bignum-len digits))
1090 (res (%allocate-bignum res-len)))
1091 (declare (type bignum-index res-len))
1092 (bignum-replace res bignum :start1 digits :end1 res-len :end2 bignum-len
1096 ;;; BIGNUM-TRUNCATE uses this to store into a bignum buffer by supplying res.
1097 ;;; When res comes in non-nil, then this foregoes allocating a result, and it
1098 ;;; normalizes the buffer instead of the would-be allocated result.
1100 ;;; We start storing into one digit higher than digits, storing a whole result
1101 ;;; digit from parts of two contiguous digits from bignum. When the loop
1102 ;;; finishes, we store the remaining bits from bignum's first digit in the
1103 ;;; first non-zero result digit, digits. We also grab some left over high
1104 ;;; bits from the last digit of bignum.
1105 (defun bignum-ashift-left-unaligned (bignum digits n-bits res-len
1106 &optional (res nil resp))
1107 (declare (type bignum-index digits res-len)
1108 (type (mod #.digit-size) n-bits))
1109 (let* ((remaining-bits (- digit-size n-bits))
1110 (res-len-1 (1- res-len))
1111 (res (or res (%allocate-bignum res-len))))
1112 (declare (type bignum-index res-len res-len-1))
1114 (j (1+ digits) (1+ j)))
1116 (setf (%bignum-ref res digits)
1117 (%ashl (%bignum-ref bignum 0) n-bits))
1118 (setf (%bignum-ref res j)
1119 (%ashr (%bignum-ref bignum i) remaining-bits))
1121 (%normalize-bignum-buffer res res-len)
1122 (%normalize-bignum res res-len)))
1123 (declare (type bignum-index i j))
1124 (setf (%bignum-ref res j)
1125 (%logior (%digit-logical-shift-right (%bignum-ref bignum i)
1127 (%ashl (%bignum-ref bignum (1+ i)) n-bits))))))
1129 ;;;; relational operators
1131 ;;; Return T iff bignum is positive.
1132 (defun bignum-plus-p (bignum)
1133 (declare (type bignum-type bignum))
1134 (%bignum-0-or-plusp bignum (%bignum-length bignum)))
1136 ;;; This compares two bignums returning -1, 0, or 1, depending on
1137 ;;; whether a is less than, equal to, or greater than b.
1138 (declaim (ftype (function (bignum bignum) (integer -1 1)) bignum-compare))
1139 (defun bignum-compare (a b)
1140 (declare (type bignum-type a b))
1141 (let* ((len-a (%bignum-length a))
1142 (len-b (%bignum-length b))
1143 (a-plusp (%bignum-0-or-plusp a len-a))
1144 (b-plusp (%bignum-0-or-plusp b len-b)))
1145 (declare (type bignum-index len-a len-b))
1146 (cond ((not (eq a-plusp b-plusp))
1149 (do ((i (1- len-a) (1- i)))
1151 (declare (type bignum-index i))
1152 (let ((a-digit (%bignum-ref a i))
1153 (b-digit (%bignum-ref b i)))
1154 (declare (type bignum-element-type a-digit b-digit))
1155 (when (%digit-greater a-digit b-digit)
1157 (when (%digit-greater b-digit a-digit)
1159 (when (zerop i) (return 0))))
1162 (t (if a-plusp -1 1)))))
1164 ;;;; float conversion
1166 ;;; Make a single or double float with the specified significand,
1167 ;;; exponent and sign.
1168 (defun single-float-from-bits (bits exp plusp)
1169 (declare (fixnum exp))
1170 (declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
1172 sb!vm:single-float-exponent-byte
1173 (logandc2 (logand #xffffffff
1174 (%bignum-ref bits 1))
1175 sb!vm:single-float-hidden-bit))))
1179 (logior res (ash -1 sb!vm:float-sign-shift))))))
1180 (defun double-float-from-bits (bits exp plusp)
1181 (declare (fixnum exp))
1182 (declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
1184 sb!vm:double-float-exponent-byte
1185 (logandc2 (ecase sb!vm::n-word-bits
1186 (32 (%bignum-ref bits 2))
1187 (64 (ash (%bignum-ref bits 1) -32)))
1188 sb!vm:double-float-hidden-bit)))
1189 (lo (logand #xffffffff (%bignum-ref bits 1))))
1190 (make-double-float (if plusp
1192 (logior hi (ash -1 sb!vm:float-sign-shift)))
1194 #!+(and long-float x86)
1195 (defun long-float-from-bits (bits exp plusp)
1196 (declare (fixnum exp))
1197 (declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
1201 (logior exp (ash 1 15)))
1202 (%bignum-ref bits 2)
1203 (%bignum-ref bits 1)))
1205 ;;; Convert Bignum to a float in the specified Format, rounding to the best
1207 (defun bignum-to-float (bignum format)
1208 (let* ((plusp (bignum-plus-p bignum))
1209 (x (if plusp bignum (negate-bignum bignum)))
1210 (len (bignum-integer-length x))
1211 (digits (float-format-digits format))
1212 (keep (+ digits digit-size))
1213 (shift (- keep len))
1214 (shifted (if (minusp shift)
1215 (bignum-ashift-right x (- shift))
1216 (bignum-ashift-left x shift)))
1217 (low (%bignum-ref shifted 0))
1218 (round-bit (ash 1 (1- digit-size))))
1219 (declare (type bignum-index len digits keep) (fixnum shift))
1220 (labels ((round-up ()
1221 (let ((rounded (add-bignums shifted round-bit)))
1222 (if (> (integer-length rounded) keep)
1223 (float-from-bits (bignum-ashift-right rounded 1)
1225 (float-from-bits rounded len))))
1226 (float-from-bits (bits len)
1227 (declare (type bignum-index len))
1230 (single-float-from-bits
1232 (check-exponent len sb!vm:single-float-bias
1233 sb!vm:single-float-normal-exponent-max)
1236 (double-float-from-bits
1238 (check-exponent len sb!vm:double-float-bias
1239 sb!vm:double-float-normal-exponent-max)
1243 (long-float-from-bits
1245 (check-exponent len sb!vm:long-float-bias
1246 sb!vm:long-float-normal-exponent-max)
1248 (check-exponent (exp bias max)
1249 (declare (type bignum-index len))
1250 (let ((exp (+ exp bias)))
1252 ;; Why a SIMPLE-TYPE-ERROR? Well, this is mainly
1253 ;; called by COERCE, which requires an error of
1254 ;; TYPE-ERROR if the conversion can't happen
1255 ;; (except in certain circumstances when we are
1256 ;; coercing to a FUNCTION) -- CSR, 2002-09-18
1257 (error 'simple-type-error
1258 :format-control "Too large to be represented as a ~S:~% ~S"
1259 :format-arguments (list format x)
1260 :expected-type format
1265 ;; Round down if round bit is 0.
1266 ((not (logtest round-bit low))
1267 (float-from-bits shifted len))
1268 ;; If only round bit is set, then round to even.
1269 ((and (= low round-bit)
1270 (dotimes (i (- (%bignum-length x) (ceiling keep digit-size))
1272 (unless (zerop (%bignum-ref x i)) (return nil))))
1273 (let ((next (%bignum-ref shifted 1)))
1276 (float-from-bits shifted len))))
1277 ;; Otherwise, round up.
1281 ;;;; integer length and logbitp/logcount
1283 (defun bignum-buffer-integer-length (bignum len)
1284 (declare (type bignum-type bignum))
1285 (let* ((len-1 (1- len))
1286 (digit (%bignum-ref bignum len-1)))
1287 (declare (type bignum-index len len-1)
1288 (type bignum-element-type digit))
1289 (+ (integer-length (%fixnum-digit-with-correct-sign digit))
1290 (* len-1 digit-size))))
1292 (defun bignum-integer-length (bignum)
1293 (declare (type bignum-type bignum))
1294 (bignum-buffer-integer-length bignum (%bignum-length bignum)))
1296 (defun bignum-logbitp (index bignum)
1297 (declare (type bignum-type bignum))
1298 (let ((len (%bignum-length bignum)))
1299 (declare (type bignum-index len))
1300 (multiple-value-bind (word-index bit-index)
1301 (floor index digit-size)
1302 (if (>= word-index len)
1303 (not (bignum-plus-p bignum))
1304 (logbitp bit-index (%bignum-ref bignum word-index))))))
1306 (defun bignum-logcount (bignum)
1307 (declare (type bignum-type bignum))
1308 (let ((length (%bignum-length bignum))
1310 (declare (type bignum-index length)
1312 (do ((index 0 (1+ index)))
1314 (if (%bignum-0-or-plusp bignum length)
1316 (- (* length digit-size) result)))
1317 (let ((digit (%bignum-ref bignum index)))
1318 (declare (type bignum-element-type digit))
1319 (incf result (logcount digit))))))
1321 ;;;; logical operations
1325 (defun bignum-logical-not (a)
1326 (declare (type bignum-type a))
1327 (let* ((len (%bignum-length a))
1328 (res (%allocate-bignum len)))
1329 (declare (type bignum-index len))
1330 (dotimes (i len res)
1331 (declare (type bignum-index i))
1332 (setf (%bignum-ref res i) (%lognot (%bignum-ref a i))))))
1336 (defun bignum-logical-and (a b)
1337 (declare (type bignum-type a b))
1338 (let* ((len-a (%bignum-length a))
1339 (len-b (%bignum-length b))
1340 (a-plusp (%bignum-0-or-plusp a len-a))
1341 (b-plusp (%bignum-0-or-plusp b len-b)))
1342 (declare (type bignum-index len-a len-b))
1346 (logand-shorter-positive a len-a b (%allocate-bignum len-a))
1347 (logand-shorter-negative a len-a b len-b (%allocate-bignum len-b))))
1350 (logand-shorter-positive b len-b a (%allocate-bignum len-b))
1351 (logand-shorter-negative b len-b a len-a (%allocate-bignum len-a))))
1352 (t (logand-shorter-positive a len-a b (%allocate-bignum len-a))))))
1354 ;;; This takes a shorter bignum, a and len-a, that is positive. Because this
1355 ;;; is AND, we don't care about any bits longer than a's since its infinite 0
1356 ;;; sign bits will mask the other bits out of b. The result is len-a big.
1357 (defun logand-shorter-positive (a len-a b res)
1358 (declare (type bignum-type a b res)
1359 (type bignum-index len-a))
1361 (declare (type bignum-index i))
1362 (setf (%bignum-ref res i)
1363 (%logand (%bignum-ref a i) (%bignum-ref b i))))
1364 (%normalize-bignum res len-a))
1366 ;;; This takes a shorter bignum, a and len-a, that is negative. Because this
1367 ;;; is AND, we just copy any bits longer than a's since its infinite 1 sign
1368 ;;; bits will include any bits from b. The result is len-b big.
1369 (defun logand-shorter-negative (a len-a b len-b res)
1370 (declare (type bignum-type a b res)
1371 (type bignum-index len-a len-b))
1373 (declare (type bignum-index i))
1374 (setf (%bignum-ref res i)
1375 (%logand (%bignum-ref a i) (%bignum-ref b i))))
1376 (do ((i len-a (1+ i)))
1378 (declare (type bignum-index i))
1379 (setf (%bignum-ref res i) (%bignum-ref b i)))
1380 (%normalize-bignum res len-b))
1384 (defun bignum-logical-ior (a b)
1385 (declare (type bignum-type a b))
1386 (let* ((len-a (%bignum-length a))
1387 (len-b (%bignum-length b))
1388 (a-plusp (%bignum-0-or-plusp a len-a))
1389 (b-plusp (%bignum-0-or-plusp b len-b)))
1390 (declare (type bignum-index len-a len-b))
1394 (logior-shorter-positive a len-a b len-b (%allocate-bignum len-b))
1395 (logior-shorter-negative a len-a b len-b (%allocate-bignum len-b))))
1398 (logior-shorter-positive b len-b a len-a (%allocate-bignum len-a))
1399 (logior-shorter-negative b len-b a len-a (%allocate-bignum len-a))))
1400 (t (logior-shorter-positive a len-a b len-b (%allocate-bignum len-a))))))
1402 ;;; This takes a shorter bignum, a and len-a, that is positive. Because this
1403 ;;; is IOR, we don't care about any bits longer than a's since its infinite
1404 ;;; 0 sign bits will mask the other bits out of b out to len-b. The result
1406 (defun logior-shorter-positive (a len-a b len-b res)
1407 (declare (type bignum-type a b res)
1408 (type bignum-index len-a len-b))
1410 (declare (type bignum-index i))
1411 (setf (%bignum-ref res i)
1412 (%logior (%bignum-ref a i) (%bignum-ref b i))))
1413 (do ((i len-a (1+ i)))
1415 (declare (type bignum-index i))
1416 (setf (%bignum-ref res i) (%bignum-ref b i)))
1417 (%normalize-bignum res len-b))
1419 ;;; This takes a shorter bignum, a and len-a, that is negative. Because this
1420 ;;; is IOR, we just copy any bits longer than a's since its infinite 1 sign
1421 ;;; bits will include any bits from b. The result is len-b long.
1422 (defun logior-shorter-negative (a len-a b len-b res)
1423 (declare (type bignum-type a b res)
1424 (type bignum-index len-a len-b))
1426 (declare (type bignum-index i))
1427 (setf (%bignum-ref res i)
1428 (%logior (%bignum-ref a i) (%bignum-ref b i))))
1429 (do ((i len-a (1+ i))
1430 (sign (%sign-digit a len-a)))
1432 (declare (type bignum-index i))
1433 (setf (%bignum-ref res i) sign))
1434 (%normalize-bignum res len-b))
1438 (defun bignum-logical-xor (a b)
1439 (declare (type bignum-type a b))
1440 (let ((len-a (%bignum-length a))
1441 (len-b (%bignum-length b)))
1442 (declare (type bignum-index len-a len-b))
1444 (bignum-logical-xor-aux a len-a b len-b (%allocate-bignum len-b))
1445 (bignum-logical-xor-aux b len-b a len-a (%allocate-bignum len-a)))))
1447 ;;; This takes the shorter of two bignums in a and len-a. Res is len-b
1448 ;;; long. Do the XOR.
1449 (defun bignum-logical-xor-aux (a len-a b len-b res)
1450 (declare (type bignum-type a b res)
1451 (type bignum-index len-a len-b))
1453 (declare (type bignum-index i))
1454 (setf (%bignum-ref res i)
1455 (%logxor (%bignum-ref a i) (%bignum-ref b i))))
1456 (do ((i len-a (1+ i))
1457 (sign (%sign-digit a len-a)))
1459 (declare (type bignum-index i))
1460 (setf (%bignum-ref res i) (%logxor sign (%bignum-ref b i))))
1461 (%normalize-bignum res len-b))
1463 ;;;; There used to be a bunch of code to implement "efficient" versions of LDB
1464 ;;;; and DPB here. But it apparently was never used, so it's been deleted.
1465 ;;;; --njf, 2007-02-04
1469 ;;; This is the original sketch of the algorithm from which I implemented this
1470 ;;; TRUNCATE, assuming both operands are bignums. I should modify this to work
1471 ;;; with the documentation on my functions, as a general introduction. I've
1472 ;;; left this here just in case someone needs it in the future. Don't look at
1473 ;;; this unless reading the functions' comments leaves you at a loss. Remember
1474 ;;; this comes from Knuth, so the book might give you the right general
1479 ;;; If X's magnitude is less than Y's, then result is 0 with remainder X.
1481 ;;; Make x and y positive, copying x if it is already positive.
1483 ;;; Shift y left until there's a 1 in the 30'th bit (most significant, non-sign
1485 ;;; Just do most sig digit to determine how much to shift whole number.
1486 ;;; Shift x this much too.
1487 ;;; Remember this initial shift count.
1489 ;;; Allocate q to be len-x minus len-y quantity plus 1.
1491 ;;; i = last digit of x.
1492 ;;; k = last digit of q.
1496 ;;; j = last digit of y.
1499 ;;; if x[i] = y[j] then g = (1- (ash 1 digit-size))
1500 ;;; else g = x[i]x[i-1]/y[j].
1503 ;;; %UNSIGNED-MULTIPLY returns b and c defined below.
1504 ;;; a = x[i-1] - (logand (* g y[j]) #xFFFFFFFF).
1505 ;;; Use %UNSIGNED-MULTIPLY taking low-order result.
1506 ;;; b = (logand (ash (* g y[j-1]) (- digit-size)) (1- (ash 1 digit-size))).
1507 ;;; c = (logand (* g y[j-1]) (1- (ash 1 digit-size))).
1509 ;;; if a > b, guess is too high
1510 ;;; g = g - 1; go back to "check guess".
1511 ;;; if a = b and c > x[i-2], guess is too high
1512 ;;; g = g - 1; go back to "check guess".
1513 ;;; GUESS IS 32-BIT NUMBER, SO USE THING TO KEEP IN SPECIAL REGISTER
1514 ;;; SAME FOR A, B, AND C.
1516 ;;; Subtract g * y from x[i - len-y+1]..x[i]. See paper for doing this in step.
1517 ;;; If x[i] < 0, guess is screwed up.
1518 ;;; negative g, then add 1
1519 ;;; zero or positive g, then subtract 1
1520 ;;; AND add y back into x[len-y+1..i].
1526 ;;; If k>=0, goto LOOP.
1528 ;;; Now quotient is good, but remainder is not.
1529 ;;; Shift x right by saved initial left shifting count.
1531 ;;; Check quotient and remainder signs.
1532 ;;; x pos y pos --> q pos r pos
1533 ;;; x pos y neg --> q neg r pos
1534 ;;; x neg y pos --> q neg r neg
1535 ;;; x neg y neg --> q pos r neg
1537 ;;; Normalize quotient and remainder. Cons result if necessary.
1540 ;;; This used to be split into multiple functions, which shared state
1541 ;;; in special variables *TRUNCATE-X* and *TRUNCATE-Y*. Having so many
1542 ;;; special variable accesses in tight inner loops was having a large
1543 ;;; effect on performance, so the helper functions have now been
1544 ;;; refactored into local functions and the special variables into
1545 ;;; lexicals. There was also a lot of boxing and unboxing of
1546 ;;; (UNSIGNED-BYTE 32)'s going on, which this refactoring
1547 ;;; eliminated. This improves the performance on some CL-BENCH tests
1548 ;;; by up to 50%, which is probably signigicant enough to justify the
1549 ;;; reduction in readability that was introduced. --JES, 2004-08-07
1550 (defun bignum-truncate (x y)
1551 (declare (type bignum-type x y))
1552 (let (truncate-x truncate-y)
1554 ;;; Divide X by Y when Y is a single bignum digit. BIGNUM-TRUNCATE
1555 ;;; fixes up the quotient and remainder with respect to sign and
1558 ;;; We don't have to worry about shifting Y to make its most
1559 ;;; significant digit sufficiently large for %FLOOR to return
1560 ;;; digit-size quantities for the q-digit and r-digit. If Y is
1561 ;;; a single digit bignum, it is already large enough for
1562 ;;; %FLOOR. That is, it has some bits on pretty high in the
1564 ((bignum-truncate-single-digit (x len-x y)
1565 (declare (type bignum-index len-x))
1566 (let ((y (%bignum-ref y 0)))
1567 (declare (type bignum-element-type y))
1568 (if (not (logtest y (1- y)))
1569 ;; Y is a power of two.
1570 ;; SHIFT-RIGHT-UNALIGNED won't do the right thing
1571 ;; with a shift count of 0 or -1, so special case this.
1573 (error 'division-by-zero))
1575 ;; We could probably get away with (VALUES X 0)
1576 ;; here, but it's not clear that some of the
1577 ;; normalization logic further down would avoid
1578 ;; mutilating X. Just go ahead and cons, consing's
1580 (values (copy-bignum x len-x) 0))
1582 (let ((n-bits (1- (integer-length y))))
1584 (shift-right-unaligned x 0 n-bits len-x
1586 (setf (%bignum-ref res j)
1587 (%ashr (%bignum-ref x i) n-bits))
1590 (logand (%bignum-ref x 0) (1- y))))))
1591 (do ((i (1- len-x) (1- i))
1592 (q (%allocate-bignum len-x))
1595 (let ((rem (%allocate-bignum 1)))
1596 (setf (%bignum-ref rem 0) r)
1598 (declare (type bignum-element-type r))
1599 (multiple-value-bind (q-digit r-digit)
1600 (%floor r (%bignum-ref x i) y)
1601 (declare (type bignum-element-type q-digit r-digit))
1602 (setf (%bignum-ref q i) q-digit)
1603 (setf r r-digit))))))
1604 ;;; This returns a guess for the next division step. Y1 is the
1605 ;;; highest y digit, and y2 is the second to highest y
1606 ;;; digit. The x... variables are the three highest x digits
1607 ;;; for the next division step.
1609 ;;; From Knuth, our guess is either all ones or x-i and x-i-1
1610 ;;; divided by y1, depending on whether x-i and y1 are the
1611 ;;; same. We test this guess by determining whether guess*y2
1612 ;;; is greater than the three high digits of x minus guess*y1
1613 ;;; shifted left one digit:
1614 ;;; ------------------------------
1615 ;;; | x-i | x-i-1 | x-i-2 |
1616 ;;; ------------------------------
1617 ;;; ------------------------------
1618 ;;; - | g*y1 high | g*y1 low | 0 |
1619 ;;; ------------------------------
1620 ;;; ... < guess*y2 ???
1621 ;;; If guess*y2 is greater, then we decrement our guess by one
1622 ;;; and try again. This returns a guess that is either
1623 ;;; correct or one too large.
1624 (bignum-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
1625 (declare (type bignum-element-type y1 y2 x-i x-i-1 x-i-2))
1626 (let ((guess (if (%digit-compare x-i y1)
1628 (%floor x-i x-i-1 y1))))
1629 (declare (type bignum-element-type guess))
1631 (multiple-value-bind (high-guess*y1 low-guess*y1)
1632 (%multiply guess y1)
1633 (declare (type bignum-element-type low-guess*y1
1635 (multiple-value-bind (high-guess*y2 low-guess*y2)
1636 (%multiply guess y2)
1637 (declare (type bignum-element-type high-guess*y2
1639 (multiple-value-bind (middle-digit borrow)
1640 (%subtract-with-borrow x-i-1 low-guess*y1 1)
1641 (declare (type bignum-element-type middle-digit)
1643 ;; Supplying borrow of 1 means there was no
1644 ;; borrow, and we know x-i-2 minus 0 requires
1646 (let ((high-digit (%subtract-with-borrow x-i
1649 (declare (type bignum-element-type high-digit))
1650 (if (and (%digit-compare high-digit 0)
1651 (or (%digit-greater high-guess*y2
1653 (and (%digit-compare middle-digit
1655 (%digit-greater low-guess*y2
1657 (setf guess (%subtract-with-borrow guess 1 1))
1658 (return guess)))))))))
1659 ;;; Divide TRUNCATE-X by TRUNCATE-Y, returning the quotient
1660 ;;; and destructively modifying TRUNCATE-X so that it holds
1663 ;;; LEN-X and LEN-Y tell us how much of the buffers we care about.
1665 ;;; TRUNCATE-X definitely has at least three digits, and it has one
1666 ;;; more than TRUNCATE-Y. This keeps i, i-1, i-2, and low-x-digit
1667 ;;; happy. Thanks to SHIFT-AND-STORE-TRUNCATE-BUFFERS.
1668 (return-quotient-leaving-remainder (len-x len-y)
1669 (declare (type bignum-index len-x len-y))
1670 (let* ((len-q (- len-x len-y))
1671 ;; Add one for extra sign digit in case high bit is on.
1672 (q (%allocate-bignum (1+ len-q)))
1674 (y1 (%bignum-ref truncate-y (1- len-y)))
1675 (y2 (%bignum-ref truncate-y (- len-y 2)))
1679 (low-x-digit (- i len-y)))
1680 (declare (type bignum-index len-q k i i-1 i-2 low-x-digit)
1681 (type bignum-element-type y1 y2))
1683 (setf (%bignum-ref q k)
1684 (try-bignum-truncate-guess
1685 ;; This modifies TRUNCATE-X. Must access
1686 ;; elements each pass.
1687 (bignum-truncate-guess y1 y2
1688 (%bignum-ref truncate-x i)
1689 (%bignum-ref truncate-x i-1)
1690 (%bignum-ref truncate-x i-2))
1692 (cond ((zerop k) (return))
1695 (shiftf i i-1 i-2 (1- i-2)))))
1697 ;;; This takes a digit guess, multiplies it by TRUNCATE-Y for a
1698 ;;; result one greater in length than LEN-Y, and subtracts this result
1699 ;;; from TRUNCATE-X. LOW-X-DIGIT is the first digit of X to start
1700 ;;; the subtraction, and we know X is long enough to subtract a LEN-Y
1701 ;;; plus one length bignum from it. Next we check the result of the
1702 ;;; subtraction, and if the high digit in X became negative, then our
1703 ;;; guess was one too big. In this case, return one less than GUESS
1704 ;;; passed in, and add one value of Y back into X to account for
1705 ;;; subtracting one too many. Knuth shows that the guess is wrong on
1706 ;;; the order of 3/b, where b is the base (2 to the digit-size power)
1707 ;;; -- pretty rarely.
1708 (try-bignum-truncate-guess (guess len-y low-x-digit)
1709 (declare (type bignum-index low-x-digit len-y)
1710 (type bignum-element-type guess))
1711 (let ((carry-digit 0)
1714 (declare (type bignum-element-type carry-digit)
1715 (type bignum-index i)
1717 ;; Multiply guess and divisor, subtracting from dividend
1720 (multiple-value-bind (high-digit low-digit)
1721 (%multiply-and-add guess
1722 (%bignum-ref truncate-y j)
1724 (declare (type bignum-element-type high-digit low-digit))
1725 (setf carry-digit high-digit)
1726 (multiple-value-bind (x temp-borrow)
1727 (%subtract-with-borrow (%bignum-ref truncate-x i)
1730 (declare (type bignum-element-type x)
1731 (fixnum temp-borrow))
1732 (setf (%bignum-ref truncate-x i) x)
1733 (setf borrow temp-borrow)))
1735 (setf (%bignum-ref truncate-x i)
1736 (%subtract-with-borrow (%bignum-ref truncate-x i)
1737 carry-digit borrow))
1738 ;; See whether guess is off by one, adding one
1739 ;; Y back in if necessary.
1740 (cond ((%digit-0-or-plusp (%bignum-ref truncate-x i))
1743 ;; If subtraction has negative result, add one
1744 ;; divisor value back in. The guess was one too
1745 ;; large in magnitude.
1746 (let ((i low-x-digit)
1749 (multiple-value-bind (v k)
1750 (%add-with-carry (%bignum-ref truncate-y j)
1751 (%bignum-ref truncate-x i)
1753 (declare (type bignum-element-type v))
1754 (setf (%bignum-ref truncate-x i) v)
1757 (setf (%bignum-ref truncate-x i)
1758 (%add-with-carry (%bignum-ref truncate-x i)
1760 (%subtract-with-borrow guess 1 1)))))
1761 ;;; This returns the amount to shift y to place a one in the
1762 ;;; second highest bit. Y must be positive. If the last digit
1763 ;;; of y is zero, then y has a one in the previous digit's
1764 ;;; sign bit, so we know it will take one less than digit-size
1765 ;;; to get a one where we want. Otherwise, we count how many
1766 ;;; right shifts it takes to get zero; subtracting this value
1767 ;;; from digit-size tells us how many high zeros there are
1768 ;;; which is one more than the shift amount sought.
1770 ;;; Note: This is exactly the same as one less than the
1771 ;;; integer-length of the last digit subtracted from the
1774 ;;; We shift y to make it sufficiently large that doing the
1775 ;;; 2*digit-size by digit-size %FLOOR calls ensures the quotient and
1776 ;;; remainder fit in digit-size.
1777 (shift-y-for-truncate (y)
1778 (let* ((len (%bignum-length y))
1779 (last (%bignum-ref y (1- len))))
1780 (declare (type bignum-index len)
1781 (type bignum-element-type last))
1782 (- digit-size (integer-length last) 1)))
1783 ;;; Stores two bignums into the truncation bignum buffers,
1784 ;;; shifting them on the way in. This assumes x and y are
1785 ;;; positive and at least two in length, and it assumes
1786 ;;; truncate-x and truncate-y are one digit longer than x and
1788 (shift-and-store-truncate-buffers (x len-x y len-y shift)
1789 (declare (type bignum-index len-x len-y)
1790 (type (integer 0 (#.digit-size)) shift))
1791 (cond ((zerop shift)
1792 (bignum-replace truncate-x x :end1 len-x)
1793 (bignum-replace truncate-y y :end1 len-y))
1795 (bignum-ashift-left-unaligned x 0 shift (1+ len-x)
1797 (bignum-ashift-left-unaligned y 0 shift (1+ len-y)
1798 truncate-y))))) ;; LABELS
1799 ;;; Divide X by Y returning the quotient and remainder. In the
1800 ;;; general case, we shift Y to set up for the algorithm, and we
1801 ;;; use two buffers to save consing intermediate values. X gets
1802 ;;; destructively modified to become the remainder, and we have
1803 ;;; to shift it to account for the initial Y shift. After we
1804 ;;; multiple bind q and r, we first fix up the signs and then
1805 ;;; return the normalized results.
1806 (let* ((x-plusp (%bignum-0-or-plusp x (%bignum-length x)))
1807 (y-plusp (%bignum-0-or-plusp y (%bignum-length y)))
1808 (x (if x-plusp x (negate-bignum x nil)))
1809 (y (if y-plusp y (negate-bignum y nil)))
1810 (len-x (%bignum-length x))
1811 (len-y (%bignum-length y)))
1812 (multiple-value-bind (q r)
1814 (bignum-truncate-single-digit x len-x y))
1815 ((plusp (bignum-compare y x))
1816 (let ((res (%allocate-bignum len-x)))
1818 (setf (%bignum-ref res i) (%bignum-ref x i)))
1821 (let ((len-x+1 (1+ len-x)))
1822 (setf truncate-x (%allocate-bignum len-x+1))
1823 (setf truncate-y (%allocate-bignum (1+ len-y)))
1824 (let ((y-shift (shift-y-for-truncate y)))
1825 (shift-and-store-truncate-buffers x len-x y
1827 (values (return-quotient-leaving-remainder len-x+1
1829 ;; Now that RETURN-QUOTIENT-LEAVING-REMAINDER
1830 ;; has executed, we just tidy up the remainder
1831 ;; (in TRUNCATE-X) and return it.
1834 (let ((res (%allocate-bignum len-y)))
1835 (declare (type bignum-type res))
1836 (bignum-replace res truncate-x :end2 len-y)
1837 (%normalize-bignum res len-y)))
1839 (shift-right-unaligned
1840 truncate-x 0 y-shift len-y
1842 (setf (%bignum-ref res j)
1843 (%ashr (%bignum-ref truncate-x i)
1845 (%normalize-bignum res res-len))
1847 (let ((quotient (cond ((eq x-plusp y-plusp) q)
1848 ((typep q 'fixnum) (the fixnum (- q)))
1849 (t (negate-bignum-in-place q))))
1850 (rem (cond (x-plusp r)
1851 ((typep r 'fixnum) (the fixnum (- r)))
1852 (t (negate-bignum-in-place r)))))
1853 (values (if (typep quotient 'fixnum)
1855 (%normalize-bignum quotient (%bignum-length quotient)))
1856 (if (typep rem 'fixnum)
1858 (%normalize-bignum rem (%bignum-length rem))))))))))
1861 ;;;; There used to be a pile of code for implementing division for bignum digits
1862 ;;;; for machines that don't have a 2*digit-size by digit-size divide instruction.
1863 ;;;; This happens to be most machines, but all the SBCL ports seem to be content
1864 ;;;; to implement SB-BIGNUM:%FLOOR as a VOP rather than using the code here.
1865 ;;;; So it's been deleted. --njf, 2007-02-04
1867 ;;;; general utilities
1869 ;;; Internal in-place operations use this to fixup remaining digits in the
1870 ;;; incoming data, such as in-place shifting. This is basically the same as
1871 ;;; the first form in %NORMALIZE-BIGNUM, but we return the length of the buffer
1872 ;;; instead of shrinking the bignum.
1873 #!-sb-fluid (declaim (sb!ext:maybe-inline %normalize-bignum-buffer))
1874 (defun %normalize-bignum-buffer (result len)
1875 (declare (type bignum-type result)
1876 (type bignum-index len))
1878 (do ((next-digit (%bignum-ref result (- len 2))
1879 (%bignum-ref result (- len 2)))
1880 (sign-digit (%bignum-ref result (1- len)) next-digit))
1881 ((not (zerop (logxor sign-digit (%ashr next-digit (1- digit-size))))))
1883 (setf (%bignum-ref result len) 0)
1888 ;;; This drops the last digit if it is unnecessary sign information. It repeats
1889 ;;; this as needed, possibly ending with a fixnum. If the resulting length from
1890 ;;; shrinking is one, see whether our one word is a fixnum. Shift the possible
1891 ;;; fixnum bits completely out of the word, and compare this with shifting the
1892 ;;; sign bit all the way through. If the bits are all 1's or 0's in both words,
1893 ;;; then there are just sign bits between the fixnum bits and the sign bit. If
1894 ;;; we do have a fixnum, shift it over for the two low-tag bits.
1895 (defun %normalize-bignum (result len)
1896 (declare (type bignum-type result)
1897 (type bignum-index len)
1898 (inline %normalize-bignum-buffer))
1899 (let ((newlen (%normalize-bignum-buffer result len)))
1900 (declare (type bignum-index newlen))
1901 (unless (= newlen len)
1902 (%bignum-set-length result newlen))
1904 (let ((digit (%bignum-ref result 0)))
1905 (if (= (%ashr digit sb!vm:n-positive-fixnum-bits)
1906 (%ashr digit (1- digit-size)))
1907 (%fixnum-digit-with-correct-sign digit)
1911 ;;; This drops the last digit if it is unnecessary sign information. It
1912 ;;; repeats this as needed, possibly ending with a fixnum magnitude but never
1913 ;;; returning a fixnum.
1914 (defun %mostly-normalize-bignum (result len)
1915 (declare (type bignum-type result)
1916 (type bignum-index len)
1917 (inline %normalize-bignum-buffer))
1918 (let ((newlen (%normalize-bignum-buffer result len)))
1919 (declare (type bignum-index newlen))
1920 (unless (= newlen len)
1921 (%bignum-set-length result newlen))
1926 ;;; the bignum case of the SXHASH function
1927 (defun sxhash-bignum (x)
1928 (let ((result 316495330))
1929 (declare (type fixnum result))
1930 (dotimes (i (%bignum-length x))
1931 (declare (type index i))
1932 (let ((xi (%bignum-ref x i)))
1934 (logand most-positive-fixnum