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
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-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 %bigfloor %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 ;;; Relational operators:
69 ;;; Note: The floating routines know about the float representation.
72 ;;; There might be a problem with various LET's and parameters that take a
73 ;;; digit value. We need to write these so those things stay in machine
74 ;;; registers and number stack slots. I bind locals to these values, and I
75 ;;; use function on them -- ZEROP, ASH, etc.
78 ;;; In shifting and byte operations, I use masks and logical operations that
79 ;;; could result in intermediate bignums. This is hidden by the current system,
80 ;;; but I may need to write these in a way that keeps these masks and logical
81 ;;; operations from diving into the Lisp level bignum code.
85 ;;; logior, logxor, logand
86 ;;; depending on relationals, < (twice) and <= (twice)
87 ;;; or write compare thing (twice).
88 ;;; LDB on fixnum with bignum result.
89 ;;; DPB on fixnum with bignum result.
90 ;;; TRUNCATE returns zero or one as one value and fixnum or minus fixnum
91 ;;; for the other value when given (truncate fixnum bignum).
92 ;;; Returns (truncate bignum fixnum) otherwise.
94 ;;; subtraction (twice)
97 ;;; Write MASK-FIELD and DEPOSIT-FIELD in terms of logical operations.
99 ;;; IF (/ x y) with bignums:
100 ;;; do the truncate, and if rem is 0, return quotient.
103 ;;; "truncate" each by gcd, ignoring remainder 0.
104 ;;; form ratio of each result, bottom is positive.
106 ;;;; What's a bignum?
108 (defconstant digit-size sb
!vm
:n-word-bits
)
110 (defconstant maximum-bignum-length
(1- (ash 1 (- sb
!vm
:n-word-bits
111 sb
!vm
:n-widetag-bits
))))
113 (defconstant all-ones-digit
(1- (ash 1 sb
!vm
:n-word-bits
)))
115 ;;;; internal inline routines
117 ;;; %ALLOCATE-BIGNUM must zero all elements.
118 (defun %allocate-bignum
(length)
119 (declare (type bignum-length length
))
120 (%allocate-bignum length
))
122 ;;; Extract the length of the bignum.
123 (defun %bignum-length
(bignum)
124 (declare (type bignum bignum
))
125 (%bignum-length bignum
))
127 ;;; %BIGNUM-REF needs to access bignums as obviously as possible, and it needs
128 ;;; to be able to return the digit somewhere no one looks for real objects.
129 (defun %bignum-ref
(bignum i
)
130 (declare (type bignum bignum
)
131 (type bignum-index i
))
132 (%bignum-ref bignum i
))
133 (defun %bignum-set
(bignum i value
)
134 (declare (type bignum bignum
)
135 (type bignum-index i
)
136 (type bignum-element-type value
))
137 (%bignum-set bignum i value
))
139 ;;; Return T if digit is positive, or NIL if negative.
140 (defun %digit-0-or-plusp
(digit)
141 (declare (type bignum-element-type digit
))
142 (not (logbitp (1- digit-size
) digit
)))
144 #!-sb-fluid
(declaim (inline %bignum-0-or-plusp
))
145 (defun %bignum-0-or-plusp
(bignum len
)
146 (declare (type bignum bignum
)
147 (type bignum-length len
))
148 (%digit-0-or-plusp
(%bignum-ref bignum
(1- len
))))
150 ;;; This should be in assembler, and should not cons intermediate
151 ;;; results. It returns a bignum digit and a carry resulting from adding
152 ;;; together a, b, and an incoming carry.
153 (defun %add-with-carry
(a b carry
)
154 (declare (type bignum-element-type a b
)
155 (type (mod 2) carry
))
156 (%add-with-carry a b carry
))
158 ;;; This should be in assembler, and should not cons intermediate
159 ;;; results. It returns a bignum digit and a borrow resulting from
160 ;;; subtracting b from a, and subtracting a possible incoming borrow.
162 ;;; We really do: a - b - 1 + borrow, where borrow is either 0 or 1.
163 (defun %subtract-with-borrow
(a b borrow
)
164 (declare (type bignum-element-type a b
)
165 (type (mod 2) borrow
))
166 (%subtract-with-borrow a b borrow
))
168 ;;; Multiply two digit-size numbers, returning a 2*digit-size result
169 ;;; split into two digit-size quantities.
170 (defun %multiply
(x y
)
171 (declare (type bignum-element-type x y
))
174 ;;; This multiplies x-digit and y-digit, producing high and low digits
175 ;;; manifesting the result. Then it adds the low digit, res-digit, and
176 ;;; carry-in-digit. Any carries (note, you still have to add two digits
177 ;;; at a time possibly producing two carries) from adding these three
178 ;;; digits get added to the high digit from the multiply, producing the
179 ;;; next carry digit. Res-digit is optional since two uses of this
180 ;;; primitive multiplies a single digit bignum by a multiple digit
181 ;;; bignum, and in this situation there is no need for a result buffer
182 ;;; accumulating partial results which is where the res-digit comes
184 (defun %multiply-and-add
(x-digit y-digit carry-in-digit
185 &optional
(res-digit 0))
186 (declare (type bignum-element-type x-digit y-digit res-digit carry-in-digit
))
187 (%multiply-and-add x-digit y-digit carry-in-digit res-digit
))
189 (defun %lognot
(digit)
190 (declare (type bignum-element-type digit
))
193 ;;; Each of these does the digit-size unsigned op.
194 (declaim (inline %logand %logior %logxor
))
196 (declare (type bignum-element-type a b
))
199 (declare (type bignum-element-type a b
))
202 (declare (type bignum-element-type a b
))
205 ;;; This takes a fixnum and sets it up as an unsigned digit-size
207 (defun %fixnum-to-digit
(x)
209 (logand x
(1- (ash 1 digit-size
))))
212 ;;; This takes three digits and returns the FLOOR'ed result of
213 ;;; dividing the first two as a 2*digit-size integer by the third.
215 ;;; Do weird LET and SETQ stuff to bamboozle the compiler into allowing
216 ;;; the %BIGFLOOR transform to expand into pseudo-assembler for which the
217 ;;; compiler can later correctly allocate registers.
218 (defun %bigfloor
(a b c
)
219 (let ((a a
) (b b
) (c c
))
220 (declare (type bignum-element-type a b c
))
224 ;;; Convert the digit to a regular integer assuming that the digit is signed.
225 (defun %fixnum-digit-with-correct-sign
(digit)
226 (declare (type bignum-element-type digit
))
227 (if (logbitp (1- digit-size
) digit
)
228 (logior digit
(ash -
1 digit-size
))
231 ;;; Do an arithmetic shift right of data even though bignum-element-type is
233 (defun %ashr
(data count
)
234 (declare (type bignum-element-type data
)
235 (type (mod #.sb
!vm
:n-word-bits
) count
))
238 ;;; This takes a digit-size quantity and shifts it to the left,
239 ;;; returning a digit-size quantity.
240 (defun %ashl
(data count
)
241 (declare (type bignum-element-type data
)
242 (type (mod #.sb
!vm
:n-word-bits
) count
))
245 ;;; Do an unsigned (logical) right shift of a digit by Count.
246 (defun %digit-logical-shift-right
(data count
)
247 (declare (type bignum-element-type data
)
248 (type (mod #.sb
!vm
:n-word-bits
) count
))
249 (%digit-logical-shift-right data count
))
251 ;;; Change the length of bignum to be newlen. Newlen must be the same or
252 ;;; smaller than the old length, and any elements beyond newlen must be zeroed.
253 (defun %bignum-set-length
(bignum newlen
)
254 (declare (type bignum bignum
)
255 (type bignum-length newlen
))
256 (%bignum-set-length bignum newlen
))
258 ;;; This returns 0 or "-1" depending on whether the bignum is positive. This
259 ;;; is suitable for infinite sign extension to complete additions,
260 ;;; subtractions, negations, etc. This cannot return a -1 represented as
261 ;;; a negative fixnum since it would then have to low zeros.
262 #!-sb-fluid
(declaim (inline %sign-digit
))
263 (defun %sign-digit
(bignum len
)
264 (declare (type bignum bignum
)
265 (type bignum-length len
))
266 (%ashr
(%bignum-ref bignum
(1- len
)) (1- digit-size
)))
268 (declaim (inline bignum-plus-p
))
269 (defun bignum-plus-p (bignum)
270 (declare (type bignum bignum
))
271 (%bignum-0-or-plusp bignum
(%bignum-length bignum
)))
273 (declaim (optimize (speed 3) (safety 0)))
277 (defun add-bignums (a b
)
278 (declare (type bignum a b
))
279 (declare (muffle-conditions compiler-note
)) ; returns lispobj, so what.
280 (let ((len-a (%bignum-length a
))
281 (len-b (%bignum-length b
)))
282 (multiple-value-bind (a len-a b len-b
)
284 (values a len-a b len-b
)
285 (values b len-b a len-a
))
286 (declare (type bignum a b
)
287 (type bignum-length len-a len-b
))
288 (let* ((len-res (1+ len-a
))
289 (res (%allocate-bignum len-res
))
291 (declare (type bignum-length len-res
)
293 (type (mod 2) carry
))
295 (declare (type bignum-index i
))
296 (multiple-value-bind (v k
)
297 (%add-with-carry
(%bignum-ref a i
) (%bignum-ref b i
) carry
)
298 (declare (type bignum-element-type v
)
300 (setf (%bignum-ref res i
) v
)
303 (finish-add a res carry
(%sign-digit b len-b
) len-b len-a
)
304 (setf (%bignum-ref res len-a
)
305 (%add-with-carry
(%sign-digit a len-a
)
306 (%sign-digit b len-b
)
308 (%normalize-bignum res len-res
)))))
310 ;;; This takes the longer of two bignums and propagates the carry through its
311 ;;; remaining high order digits.
312 (defun finish-add (a res carry sign-digit-b start end
)
313 (declare (type bignum a res
)
315 (type bignum-element-type sign-digit-b
)
316 (type bignum-index start
)
317 (type bignum-length end
))
318 (do ((i start
(1+ i
)))
320 (setf (%bignum-ref res end
)
321 (%add-with-carry
(%sign-digit a end
) sign-digit-b carry
)))
322 (declare (type bignum-index i
))
323 (multiple-value-bind (v k
)
324 (%add-with-carry
(%bignum-ref a i
) sign-digit-b carry
)
325 (setf (%bignum-ref res i
) v
)
331 (eval-when (:compile-toplevel
:execute
)
333 ;;; This subtracts b from a plugging result into res. Return-fun is the
334 ;;; function to call that fixes up the result returning any useful values, such
335 ;;; as the result. This macro may evaluate its arguments more than once.
336 (sb!xc
:defmacro subtract-bignum-loop
(a len-a b len-b res len-res return-fun
)
337 (with-unique-names (borrow a-digit a-sign b-digit b-sign i v k
)
339 (,a-sign
(%sign-digit
,a
,len-a
))
340 (,b-sign
(%sign-digit
,b
,len-b
)))
341 (declare (type bignum-element-type
,a-sign
,b-sign
))
342 (dotimes (,i
,len-res
)
343 (declare (type bignum-index
,i
))
344 (let ((,a-digit
(if (< ,i
,len-a
) (%bignum-ref
,a
,i
) ,a-sign
))
345 (,b-digit
(if (< ,i
,len-b
) (%bignum-ref
,b
,i
) ,b-sign
)))
346 (declare (type bignum-element-type
,a-digit
,b-digit
))
347 (multiple-value-bind (,v
,k
)
348 (%subtract-with-borrow
,a-digit
,b-digit
,borrow
)
349 (setf (%bignum-ref
,res
,i
) ,v
)
351 (,return-fun
,res
,len-res
))))
355 (defun subtract-bignum (a b
)
356 (declare (type bignum a b
))
357 (let* ((len-a (%bignum-length a
))
358 (len-b (%bignum-length b
))
359 (len-res (1+ (max len-a len-b
)))
360 (res (%allocate-bignum len-res
)))
361 (declare (type bignum-length len-a len-b len-res
)) ;Test len-res for bounds?
362 (subtract-bignum-loop a len-a b len-b res len-res %normalize-bignum
)))
364 ;;; Operations requiring a subtraction without the overhead of intermediate
365 ;;; results, such as GCD, use this. It assumes Result is big enough for the
367 (defun subtract-bignum-buffers-with-len (a len-a b len-b result len-res
)
368 (declare (type bignum a b result
)
369 (type bignum-length len-a len-b len-res
))
370 (subtract-bignum-loop a len-a b len-b result len-res
371 %normalize-bignum-buffer
))
373 (defun subtract-bignum-buffers (a len-a b len-b result
)
374 (declare (type bignum a b result
)
375 (type bignum-length len-a len-b
))
376 (subtract-bignum-loop a len-a b len-b result
(max len-a len-b
)
377 %normalize-bignum-buffer
))
381 (defun multiply-bignums (a b
)
382 (declare (type bignum a b
))
383 (let* ((a-plusp (bignum-plus-p a
))
384 (b-plusp (bignum-plus-p b
))
385 (a (if a-plusp a
(negate-bignum a
)))
386 (b (if b-plusp b
(negate-bignum b
)))
387 (len-a (%bignum-length a
))
388 (len-b (%bignum-length b
))
389 (len-res (+ len-a len-b
))
390 (res (%allocate-bignum len-res
))
391 (negate-res (not (eq a-plusp b-plusp
))))
392 (declare (type bignum-length len-a len-b len-res
))
394 (declare (type bignum-index i
))
395 (let ((carry-digit 0)
396 (x (%bignum-ref a i
))
398 (declare (type bignum-index k
)
399 (type bignum-element-type carry-digit x
))
401 (multiple-value-bind (big-carry res-digit
)
406 (declare (type bignum-element-type big-carry res-digit
))
407 (setf (%bignum-ref res k
) res-digit
)
408 (setf carry-digit big-carry
)
410 (setf (%bignum-ref res k
) carry-digit
)))
411 (when negate-res
(negate-bignum-in-place res
))
412 (%normalize-bignum res len-res
)))
414 (defun multiply-bignum-and-fixnum (bignum fixnum
)
415 (declare (type bignum bignum
) (type fixnum fixnum
))
416 (let* ((bignum-plus-p (bignum-plus-p bignum
))
417 (fixnum-plus-p (not (minusp fixnum
)))
418 (bignum (if bignum-plus-p bignum
(negate-bignum bignum
)))
419 (bignum-len (%bignum-length bignum
))
420 (fixnum (if fixnum-plus-p fixnum
(- fixnum
)))
421 (result (%allocate-bignum
(1+ bignum-len
)))
423 (declare (type bignum bignum result
)
424 (type bignum-element-type fixnum carry-digit
))
425 (dotimes (index bignum-len
)
426 (declare (type bignum-index index
))
427 (multiple-value-bind (next-digit low
)
428 (%multiply-and-add
(%bignum-ref bignum index
) fixnum carry-digit
)
429 (declare (type bignum-element-type next-digit low
))
430 (setf carry-digit next-digit
)
431 (setf (%bignum-ref result index
) low
)))
432 (setf (%bignum-ref result bignum-len
) carry-digit
)
433 (unless (eq bignum-plus-p fixnum-plus-p
)
434 (negate-bignum-in-place result
))
435 (%normalize-bignum result
(1+ bignum-len
))))
437 (defun multiply-fixnums (a b
)
438 (declare (fixnum a b
))
439 (declare (muffle-conditions compiler-note
)) ; returns lispobj, so what.
440 (let* ((a-minusp (minusp a
))
441 (b-minusp (minusp b
)))
442 (multiple-value-bind (high low
)
443 (%multiply
(if a-minusp
(- a
) a
)
444 (if b-minusp
(- b
) b
))
445 (declare (type bignum-element-type high low
))
446 (if (and (zerop high
)
447 (%digit-0-or-plusp low
))
448 (let ((low (truly-the (unsigned-byte #.
(1- sb
!vm
:n-word-bits
))
449 (%fixnum-digit-with-correct-sign low
))))
450 (if (eq a-minusp b-minusp
)
453 (let ((res (%allocate-bignum
2)))
454 (%bignum-set res
0 low
)
455 (%bignum-set res
1 high
)
456 (unless (eq a-minusp b-minusp
) (negate-bignum-in-place res
))
457 (%normalize-bignum res
2))))))
459 ;;;; BIGNUM-REPLACE and WITH-BIGNUM-BUFFERS
461 (eval-when (:compile-toplevel
:execute
)
463 (sb!xc
:defmacro bignum-replace
(dest
471 (sb!int
:once-only
((n-dest dest
)
473 (with-unique-names (n-start1 n-end1 n-start2 n-end2 i1 i2
)
474 (let ((end1 (or end1
`(%bignum-length
,n-dest
)))
475 (end2 (or end2
`(%bignum-length
,n-src
))))
477 `(let ((,n-start1
,start1
)
479 (do ((,i1
(1- ,end1
) (1- ,i1
))
480 (,i2
(1- ,end2
) (1- ,i2
)))
481 ((or (< ,i1
,n-start1
) (< ,i2
,n-start2
)))
482 (declare (fixnum ,i1
,i2
))
483 (%bignum-set
,n-dest
,i1
(%bignum-ref
,n-src
,i2
))))
484 (if (eql start1 start2
)
485 `(let ((,n-end1
(min ,end1
,end2
)))
486 (do ((,i1
,start1
(1+ ,i1
)))
488 (declare (type bignum-index
,i1
))
489 (%bignum-set
,n-dest
,i1
(%bignum-ref
,n-src
,i1
))))
490 `(let ((,n-end1
,end1
)
492 (do ((,i1
,start1
(1+ ,i1
))
493 (,i2
,start2
(1+ ,i2
)))
494 ((or (>= ,i1
,n-end1
) (>= ,i2
,n-end2
)))
495 (declare (type bignum-index
,i1
,i2
))
496 (%bignum-set
,n-dest
,i1
(%bignum-ref
,n-src
,i2
))))))))))
498 (sb!xc
:defmacro with-bignum-buffers
(specs &body body
)
499 "WITH-BIGNUM-BUFFERS ({(var size [init])}*) Form*"
500 (sb!int
:collect
((binds)
503 (let ((name (first spec
))
504 (size (second spec
)))
505 (binds `(,name
(%allocate-bignum
,size
)))
506 (let ((init (third spec
)))
508 (inits `(bignum-replace ,name
,init
))))))
517 (eval-when (:compile-toplevel
:load-toplevel
:execute
)
518 ;; The asserts in the GCD implementation are way too expensive to
519 ;; check in normal use, and are disabled here.
520 (sb!xc
:defmacro gcd-assert
(&rest args
)
521 (declare (ignore args
))
522 #+sb-bignum-assertions
`(assert ,@args
))
523 ;; We'll be doing a lot of modular arithmetic.
524 (sb!xc
:defmacro modularly
(form)
525 `(logand all-ones-digit
,form
)))
527 ;;; I'm not sure why I need this FTYPE declaration. Compiled by the
528 ;;; target compiler, it can deduce the return type fine, but without
529 ;;; it, we pay a heavy price in BIGNUM-GCD when compiled by the
530 ;;; cross-compiler. -- CSR, 2004-07-19
531 (declaim (ftype (sfunction (bignum bignum-length bignum bignum-length
)
532 (and unsigned-byte fixnum
))
533 bignum-factors-of-two
))
534 (defun bignum-factors-of-two (a len-a b len-b
)
535 (declare (type bignum-length len-a len-b
) (type bignum a b
))
537 (end (min len-a len-b
)))
538 ((= i end
) (error "Unexpected zero bignums?"))
539 (declare (type bignum-index i
)
540 (type bignum-length end
))
541 (let ((or-digits (%logior
(%bignum-ref a i
) (%bignum-ref b i
))))
542 (unless (zerop or-digits
)
543 (return (do ((j 0 (1+ j
))
544 (or-digits or-digits
(%ashr or-digits
1)))
545 ((oddp or-digits
) (+ (* i digit-size
) j
))
546 (declare (type (mod #.sb
!vm
:n-word-bits
) j
))))))))
548 ;;; Multiply a bignum buffer with a fixnum or a digit, storing the
549 ;;; result in another bignum buffer, and without using any
550 ;;; temporaries. Inlined to avoid boxing smallnum if it's actually a
551 ;;; digit. Needed by GCD, should possibly OAOO with
552 ;;; MULTIPLY-BIGNUM-AND-FIXNUM.
553 (declaim (inline multiply-bignum-buffer-and-smallnum-to-buffer
))
554 (defun multiply-bignum-buffer-and-smallnum-to-buffer (bignum bignum-len
556 (declare (type bignum bignum
))
557 (let* ((bignum-plus-p (%bignum-0-or-plusp bignum bignum-len
))
558 (smallnum-plus-p (not (minusp smallnum
)))
559 (smallnum (if smallnum-plus-p smallnum
(- smallnum
)))
561 (declare (type bignum bignum res
)
562 (type bignum-length bignum-len
)
563 (type bignum-element-type smallnum carry-digit
))
564 (unless bignum-plus-p
565 (negate-bignum-buffer-in-place bignum bignum-len
))
566 (dotimes (index bignum-len
)
567 (declare (type bignum-index index
))
568 (multiple-value-bind (next-digit low
)
569 (%multiply-and-add
(%bignum-ref bignum index
)
572 (declare (type bignum-element-type next-digit low
))
573 (setf carry-digit next-digit
)
574 (setf (%bignum-ref res index
) low
)))
575 (setf (%bignum-ref res bignum-len
) carry-digit
)
576 (unless bignum-plus-p
577 (negate-bignum-buffer-in-place bignum bignum-len
))
578 (let ((res-len (%normalize-bignum-buffer res
(1+ bignum-len
))))
579 (unless (eq bignum-plus-p smallnum-plus-p
)
580 (negate-bignum-buffer-in-place res res-len
))
583 ;;; Given U and V, return U / V mod 2^32. Implements the algorithm in the
584 ;;; paper, but uses some clever bit-twiddling nicked from Nickle to do it.
585 (declaim (inline bmod
))
587 (declare (muffle-conditions compiler-note
)) ; returns lispobj, so what.
588 (let ((ud (%bignum-ref u
0))
589 (vd (%bignum-ref v
0))
593 (declare (type (unsigned-byte #.sb
!vm
:n-word-bits
) ud vd umask imask m
))
594 (dotimes (i digit-size
)
595 (setf umask
(logior umask imask
))
596 (when (logtest ud umask
)
597 (setf ud
(modularly (- ud vd
)))
598 (setf m
(modularly (logior m imask
))))
599 (setf imask
(modularly (ash imask
1)))
600 (setf vd
(modularly (ash vd
1))))
603 (defun dmod (u u-len v v-len tmp1
)
604 (loop while
(> (bignum-buffer-integer-length u u-len
)
605 (+ (bignum-buffer-integer-length v v-len
)
608 (unless (zerop (%bignum-ref u
0))
609 (let* ((bmod (bmod u v
))
610 (tmp1-len (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
613 (setf u-len
(subtract-bignum-buffers u u-len
616 (bignum-abs-buffer u u-len
)))
617 (gcd-assert (zerop (%bignum-ref u
0)))
618 (setf u-len
(bignum-buffer-ashift-right u u-len digit-size
)))
619 (let* ((d (+ 1 (- (bignum-buffer-integer-length u u-len
)
620 (bignum-buffer-integer-length v v-len
))))
622 (declare (type (unsigned-byte #.
(integer-length #.sb
!vm
:n-word-bits
)) d
)
623 (type (unsigned-byte #.sb
!vm
:n-word-bits
) n
))
624 (gcd-assert (>= d
0))
625 (when (logtest (%bignum-ref u
0) n
)
627 (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
631 (setf u-len
(subtract-bignum-buffers u u-len
634 (bignum-abs-buffer u u-len
)))
637 (defconstant lower-ones-digit
(1- (ash 1 (truncate sb
!vm
:n-word-bits
2))))
639 ;;; Find D and N such that (LOGAND ALL-ONES-DIGIT (- (* D X) (* N Y))) is 0,
640 ;;; (< 0 N LOWER-ONES-DIGIT) and (< 0 (ABS D) LOWER-ONES-DIGIT).
641 (defun reduced-ratio-mod (x y
)
642 (let* ((c (bmod x y
))
645 (n2 (modularly (1+ (modularly (lognot n1
)))))
647 (declare (type (unsigned-byte #.sb
!vm
:n-word-bits
) n1 d1 n2 d2
))
648 (loop while
(> n2
(expt 2 (truncate digit-size
2))) do
649 (loop for i of-type
(mod #.sb
!vm
:n-word-bits
)
650 downfrom
(- (integer-length n1
) (integer-length n2
))
652 (when (>= n1
(modularly (ash n2 i
)))
653 (psetf n1
(modularly (- n1
(modularly (ash n2 i
))))
654 d1
(modularly (- d1
(modularly (ash d2 i
)))))))
659 (values n2
(if (>= d2
(expt 2 (1- digit-size
)))
660 (lognot (logand most-positive-fixnum
(lognot d2
)))
661 (logand lower-ones-digit d2
)))))
664 (defun copy-bignum (a &optional
(len (%bignum-length a
)))
665 (let ((b (%allocate-bignum len
)))
667 (%bignum-set-length b len
)
670 ;;; Allocate a single word bignum that holds fixnum. This is useful when
671 ;;; we are trying to mix fixnum and bignum operands.
672 #!-sb-fluid
(declaim (inline make-small-bignum
))
673 (defun make-small-bignum (fixnum)
674 (let ((res (%allocate-bignum
1)))
675 (setf (%bignum-ref res
0) (%fixnum-to-digit fixnum
))
678 ;; When the larger number is less than this many bignum digits long, revert
680 (defparameter *accelerated-gcd-cutoff
* 3)
682 ;;; Alternate between k-ary reduction with the help of
683 ;;; REDUCED-RATIO-MOD and digit modulus reduction via DMOD. Once the
684 ;;; arguments get small enough, drop through to BIGNUM-MOD-GCD (since
685 ;;; k-ary reduction can introduce spurious factors, which need to be
686 ;;; filtered out). Reference: Kenneth Weber, "The accelerated integer
687 ;;; GCD algorithm", ACM Transactions on Mathematical Software, volume
688 ;;; 21, number 1, March 1995, epp. 111-122.
689 (defun bignum-gcd (u0 v0
)
690 (declare (type bignum u0 v0
))
691 (let* ((u1 (if (bignum-plus-p u0
)
693 (negate-bignum u0 nil
)))
694 (v1 (if (bignum-plus-p v0
)
696 (negate-bignum v0 nil
))))
698 (return-from bignum-gcd u1
))
701 (let ((n (mod v1 u1
)))
702 (setf v1
(if (fixnump n
)
703 (make-small-bignum n
)
705 (if (and (= 1 (%bignum-length v1
))
706 (zerop (%bignum-ref v1
0)))
707 (return-from bignum-gcd
(%normalize-bignum u1
708 (%bignum-length u1
))))
709 (let* ((buffer-len (+ 2 (%bignum-length u1
)))
710 (u (%allocate-bignum buffer-len
))
711 (u-len (%bignum-length u1
))
712 (v (%allocate-bignum buffer-len
))
713 (v-len (%bignum-length v1
))
714 (tmp1 (%allocate-bignum buffer-len
))
716 (tmp2 (%allocate-bignum buffer-len
))
719 (bignum-factors-of-two u1
(%bignum-length u1
)
720 v1
(%bignum-length v1
))))
721 (declare (type (or null bignum-length
)
722 buffer-len u-len v-len tmp1-len tmp2-len
))
723 (bignum-replace u u1
)
724 (bignum-replace v v1
)
726 (make-gcd-bignum-odd u
727 (bignum-buffer-ashift-right u u-len
730 (make-gcd-bignum-odd v
731 (bignum-buffer-ashift-right v v-len
733 (loop until
(or (< u-len
*accelerated-gcd-cutoff
*)
737 (zerop (%bignum-ref v
0))))
739 (gcd-assert (= buffer-len
(%bignum-length u
)
741 (%bignum-length tmp1
)
742 (%bignum-length tmp2
)))
743 (if (> (bignum-buffer-integer-length u u-len
)
744 (+ #.
(truncate sb
!vm
:n-word-bits
4)
745 (bignum-buffer-integer-length v v-len
)))
746 (setf u-len
(dmod u u-len
749 (multiple-value-bind (n d
) (reduced-ratio-mod u v
)
751 (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
754 (multiply-bignum-buffer-and-smallnum-to-buffer u u-len
756 (gcd-assert (= (copy-bignum tmp2 tmp2-len
)
757 (* (copy-bignum u u-len
) d
)))
758 (gcd-assert (= (copy-bignum tmp1 tmp1-len
)
759 (* (copy-bignum v v-len
) n
)))
761 (subtract-bignum-buffers-with-len tmp1 tmp1-len
766 (gcd-assert (or (zerop (- (copy-bignum tmp1 tmp1-len
)
767 (copy-bignum tmp2 tmp2-len
)))
768 (= (copy-bignum u u-len
)
769 (- (copy-bignum tmp1 tmp1-len
)
770 (copy-bignum tmp2 tmp2-len
)))))
771 (bignum-abs-buffer u u-len
)
772 (gcd-assert (zerop (modularly u
)))))
773 (setf u-len
(make-gcd-bignum-odd u u-len
))
775 (rotatef u-len v-len
))
776 (bignum-abs-buffer u u-len
)
777 (setf u
(copy-bignum u u-len
))
778 (let ((n (bignum-mod-gcd v1 u
)))
779 (ash (bignum-mod-gcd u1
(if (fixnump n
)
780 (make-small-bignum n
)
784 (defun bignum-mod-gcd (a b
)
785 (declare (type bignum a b
))
788 ;; While the length difference of A and B is sufficiently large,
789 ;; reduce using MOD (slowish, but it should equalize the sizes of
790 ;; A and B pretty quickly). After that, use the binary GCD
791 ;; algorithm to handle the rest.
792 (loop until
(and (= (%bignum-length b
) 1) (zerop (%bignum-ref b
0))) do
793 (when (<= (%bignum-length a
) (1+ (%bignum-length b
)))
794 (return-from bignum-mod-gcd
(bignum-binary-gcd a b
)))
795 (let ((rem (mod a b
)))
797 (setf a
(make-small-bignum rem
))
800 (if (= (%bignum-length a
) 1)
801 (%normalize-bignum a
1)
804 (defun bignum-binary-gcd (a b
)
805 (declare (type bignum a b
))
806 (let* ((len-a (%bignum-length a
))
807 (len-b (%bignum-length b
)))
808 (with-bignum-buffers ((a-buffer len-a a
)
810 (res-buffer (max len-a len-b
)))
811 (let* ((factors-of-two
812 (bignum-factors-of-two a-buffer len-a
814 (len-a (make-gcd-bignum-odd
816 (bignum-buffer-ashift-right a-buffer len-a
818 (len-b (make-gcd-bignum-odd
820 (bignum-buffer-ashift-right b-buffer len-b
822 (declare (type bignum-length len-a len-b
))
829 (multiple-value-bind (u v len-v r len-r
)
830 (bignum-gcd-order-and-subtract x len-x y len-y z
)
831 (declare (type bignum-length len-v len-r
))
832 (when (and (= len-r
1) (zerop (%bignum-ref r
0)))
833 (if (zerop factors-of-two
)
834 (let ((ret (%allocate-bignum len-v
)))
836 (setf (%bignum-ref ret i
) (%bignum-ref v i
)))
837 (return (%normalize-bignum ret len-v
)))
838 (return (bignum-ashift-left v factors-of-two len-v
))))
839 (setf x v len-x len-v
)
840 (setf y r len-y
(make-gcd-bignum-odd r len-r
))
843 (defun bignum-gcd-order-and-subtract (a len-a b len-b res
)
844 (declare (type bignum-length len-a len-b
) (type bignum a b
))
845 (cond ((= len-a len-b
)
846 (do ((i (1- len-a
) (1- i
)))
848 (setf (%bignum-ref res
0) 0)
849 (values a b len-b res
1))
850 (let ((a-digit (%bignum-ref a i
))
851 (b-digit (%bignum-ref b i
)))
852 (cond ((= a-digit b-digit
))
855 (values a b len-b res
856 (subtract-bignum-buffers a len-a b len-b
860 (values b a len-a res
861 (subtract-bignum-buffers b len-b
865 (values a b len-b res
866 (subtract-bignum-buffers a len-a b len-b res
)))
868 (values b a len-a res
869 (subtract-bignum-buffers b len-b a len-a res
)))))
871 (defun make-gcd-bignum-odd (a len-a
)
872 (declare (type bignum a
) (type bignum-length len-a
))
873 (dotimes (index len-a
)
874 (declare (type bignum-index index
))
875 (do ((digit (%bignum-ref a index
) (%ashr digit
1))
876 (increment 0 (1+ increment
)))
878 (declare (type (mod #.sb
!vm
:n-word-bits
) increment
))
880 (return-from make-gcd-bignum-odd
881 (bignum-buffer-ashift-right a len-a
882 (+ (* index digit-size
)
888 (eval-when (:compile-toplevel
:execute
)
890 ;;; This negates bignum-len digits of bignum, storing the resulting digits into
891 ;;; result (possibly EQ to bignum) and returning whatever end-carry there is.
892 (sb!xc
:defmacro bignum-negate-loop
893 (bignum bignum-len
&optional
(result nil resultp
))
894 (with-unique-names (carry end value last
)
895 `(let* (,@(if (not resultp
) `(,last
))
897 (multiple-value-bind (,value
,carry
)
898 (%add-with-carry
(%lognot
(%bignum-ref
,bignum
0)) 1 0)
900 `(setf (%bignum-ref
,result
0) ,value
)
901 `(setf ,last
,value
))
905 (declare (type bit
,carry
)
906 (type bignum-index i
)
907 (type bignum-length
,end
))
909 (when (= i
,end
) (return))
910 (multiple-value-bind (,value temp
)
911 (%add-with-carry
(%lognot
(%bignum-ref
,bignum i
)) 0 ,carry
)
913 `(setf (%bignum-ref
,result i
) ,value
)
914 `(setf ,last
,value
))
917 ,(if resultp carry
`(values ,carry
,last
)))))
921 ;;; Fully-normalize is an internal optional. It cause this to always return
922 ;;; a bignum, without any extraneous digits, and it never returns a fixnum.
923 (defun negate-bignum (x &optional
(fully-normalize t
))
924 (declare (type bignum x
))
925 (let* ((len-x (%bignum-length x
))
927 (res (%allocate-bignum len-res
)))
928 (declare (type bignum-length len-x len-res
)) ;Test len-res for range?
929 (let ((carry (bignum-negate-loop x len-x res
)))
930 (setf (%bignum-ref res len-x
)
931 (%add-with-carry
(%lognot
(%sign-digit x len-x
)) 0 carry
)))
933 (%normalize-bignum res len-res
)
934 (%mostly-normalize-bignum res len-res
))))
936 ;;; This assumes bignum is positive; that is, the result of negating it will
937 ;;; stay in the provided allocated bignum.
938 (declaim (maybe-inline negate-bignum-buffer-in-place
))
939 (defun negate-bignum-buffer-in-place (bignum bignum-len
)
940 (bignum-negate-loop bignum bignum-len bignum
)
943 (defun negate-bignum-in-place (bignum)
944 (declare (inline negate-bignum-buffer-in-place
))
945 (negate-bignum-buffer-in-place bignum
(%bignum-length bignum
)))
947 (defun bignum-abs-buffer (bignum len
)
948 (unless (%bignum-0-or-plusp bignum len
)
949 (negate-bignum-buffer-in-place bignum len
)))
953 (eval-when (:compile-toplevel
:execute
)
955 ;;; This macro is used by BIGNUM-ASHIFT-RIGHT, BIGNUM-BUFFER-ASHIFT-RIGHT, and
956 ;;; BIGNUM-LDB-BIGNUM-RES. They supply a termination form that references
957 ;;; locals established by this form. Source is the source bignum. Start-digit
958 ;;; is the first digit in source from which we pull bits. Start-pos is the
959 ;;; first bit we want. Res-len-form is the form that computes the length of
960 ;;; the resulting bignum. Termination is a DO termination form with a test and
961 ;;; body. When result is supplied, it is the variable to which this binds a
962 ;;; newly allocated bignum.
964 ;;; Given start-pos, 1-31 inclusively, of shift, we form the j'th resulting
965 ;;; digit from high bits of the i'th source digit and the start-pos number of
966 ;;; bits from the i+1'th source digit.
967 (sb!xc
:defmacro shift-right-unaligned
(source
973 `(let* ((high-bits-in-first-digit (- digit-size
,start-pos
))
974 (res-len ,res-len-form
)
975 (res-len-1 (1- res-len
))
976 ,@(if result
`((,result
(%allocate-bignum res-len
)))))
977 (declare (type bignum-length res-len res-len-1
))
978 (do ((i ,start-digit
(1+ i
))
981 (declare (type bignum-index i j
))
982 (setf (%bignum-ref
,(if result result source
) j
)
983 (%logior
(%digit-logical-shift-right
(%bignum-ref
,source i
)
985 (%ashl
(%bignum-ref
,source
(1+ i
))
986 high-bits-in-first-digit
))))))
990 ;;; First compute the number of whole digits to shift, shifting them by
991 ;;; skipping them when we start to pick up bits, and the number of bits to
992 ;;; shift the remaining digits into place. If the number of digits is greater
993 ;;; than the length of the bignum, then the result is either 0 or -1. If we
994 ;;; shift on a digit boundary (that is, n-bits is zero), then we just copy
995 ;;; digits. The last branch handles the general case which uses a macro that a
996 ;;; couple other routines use. The fifth argument to the macro references
997 ;;; locals established by the macro.
998 (defun bignum-ashift-right (bignum count
)
999 (declare (type bignum bignum
)
1000 (type unsigned-byte count
))
1001 (let ((bignum-len (%bignum-length bignum
)))
1002 (cond ((fixnump count
)
1003 (multiple-value-bind (digits n-bits
) (truncate count digit-size
)
1004 (declare (type bignum-length digits
))
1006 ((>= digits bignum-len
)
1007 (if (%bignum-0-or-plusp bignum bignum-len
) 0 -
1))
1009 (bignum-ashift-right-digits bignum digits
))
1011 (shift-right-unaligned bignum digits n-bits
(- bignum-len digits
)
1013 (setf (%bignum-ref res j
)
1014 (%ashr
(%bignum-ref bignum i
) n-bits
))
1015 (%normalize-bignum res res-len
))
1017 ((> count bignum-len
)
1018 (if (%bignum-0-or-plusp bignum bignum-len
) 0 -
1))
1019 ;; Since a FIXNUM should be big enough to address anything in
1020 ;; memory, including arrays of bits, and since arrays of bits
1021 ;; take up about the same space as corresponding fixnums, there
1022 ;; should be no way that we fall through to this case: any shift
1023 ;; right by a bignum should give zero. But let's check anyway:
1024 (t (error "bignum overflow: can't shift right by ~S" count
)))))
1026 (defun bignum-ashift-right-digits (bignum digits
)
1027 (declare (type bignum bignum
)
1028 (type bignum-length digits
))
1029 (let* ((res-len (- (%bignum-length bignum
) digits
))
1030 (res (%allocate-bignum res-len
)))
1031 (declare (type bignum-length res-len
)
1033 (bignum-replace res bignum
:start2 digits
)
1034 (%normalize-bignum res res-len
)))
1036 ;;; GCD uses this for an in-place shifting operation. This is different enough
1037 ;;; from BIGNUM-ASHIFT-RIGHT that it isn't worth folding the bodies into a
1038 ;;; macro, but they share the basic algorithm. This routine foregoes a first
1039 ;;; test for digits being greater than or equal to bignum-len since that will
1040 ;;; never happen for its uses in GCD. We did fold the last branch into a macro
1041 ;;; since it was duplicated a few times, and the fifth argument to it
1042 ;;; references locals established by the macro.
1043 (defun bignum-buffer-ashift-right (bignum bignum-len x
)
1044 (declare (type bignum-length bignum-len
) (fixnum x
))
1045 (multiple-value-bind (digits n-bits
) (truncate x digit-size
)
1046 (declare (type bignum-length digits
))
1049 (let ((new-end (- bignum-len digits
)))
1050 (bignum-replace bignum bignum
:end1 new-end
:start2 digits
1052 (%normalize-bignum-buffer bignum new-end
)))
1054 (shift-right-unaligned bignum digits n-bits
(- bignum-len digits
)
1056 (setf (%bignum-ref bignum j
)
1057 (%ashr
(%bignum-ref bignum i
) n-bits
))
1058 (%normalize-bignum-buffer bignum res-len
)))))))
1060 ;;; This handles shifting a bignum buffer to provide fresh bignum data for some
1061 ;;; internal routines. We know bignum is safe when called with bignum-len.
1062 ;;; First we compute the number of whole digits to shift, shifting them
1063 ;;; starting to store farther along the result bignum. If we shift on a digit
1064 ;;; boundary (that is, n-bits is zero), then we just copy digits. The last
1065 ;;; branch handles the general case.
1066 (defun bignum-ashift-left (bignum x
&optional bignum-len
)
1067 (declare (type bignum bignum
)
1068 (type unsigned-byte x
)
1069 (type (or null bignum-length
) bignum-len
))
1071 (multiple-value-bind (digits n-bits
) (truncate x digit-size
)
1072 (let* ((bignum-len (or bignum-len
(%bignum-length bignum
)))
1073 (res-len (+ digits bignum-len
1)))
1074 (when (> res-len maximum-bignum-length
)
1075 (error "can't represent result of left shift"))
1077 (bignum-ashift-left-digits bignum bignum-len digits
)
1078 (bignum-ashift-left-unaligned bignum digits n-bits res-len
))))
1079 ;; Left shift by a number too big to be represented as a fixnum
1080 ;; would exceed our memory capacity, since a fixnum is big enough
1081 ;; to index any array, including a bit array.
1082 (error "can't represent result of left shift")))
1084 (defun bignum-ashift-left-digits (bignum bignum-len digits
)
1085 (declare (type bignum-length bignum-len digits
))
1086 (let* ((res-len (+ bignum-len digits
))
1087 (res (%allocate-bignum res-len
)))
1088 (declare (type bignum-length res-len
))
1089 (bignum-replace res bignum
:start1 digits
:end1 res-len
:end2 bignum-len
1093 ;;; BIGNUM-TRUNCATE uses this to store into a bignum buffer by supplying res.
1094 ;;; When res comes in non-nil, then this foregoes allocating a result, and it
1095 ;;; normalizes the buffer instead of the would-be allocated result.
1097 ;;; We start storing into one digit higher than digits, storing a whole result
1098 ;;; digit from parts of two contiguous digits from bignum. When the loop
1099 ;;; finishes, we store the remaining bits from bignum's first digit in the
1100 ;;; first non-zero result digit, digits. We also grab some left over high
1101 ;;; bits from the last digit of bignum.
1102 (defun bignum-ashift-left-unaligned (bignum digits n-bits res-len
1103 &optional
(res nil resp
))
1104 (declare (type bignum-length digits res-len
)
1105 (type (mod #.digit-size
) n-bits
))
1106 (let* ((remaining-bits (- digit-size n-bits
))
1107 (res-len-1 (1- res-len
))
1108 (res (or res
(%allocate-bignum res-len
))))
1109 (declare (type bignum-length res-len res-len-1
))
1111 (j (1+ digits
) (1+ j
)))
1113 (setf (%bignum-ref res digits
)
1114 (%ashl
(%bignum-ref bignum
0) n-bits
))
1115 (setf (%bignum-ref res j
)
1116 (%ashr
(%bignum-ref bignum i
) remaining-bits
))
1118 (%normalize-bignum-buffer res res-len
)
1119 (%normalize-bignum res res-len
)))
1120 (declare (type bignum-index i j
))
1121 (setf (%bignum-ref res j
)
1122 (%logior
(%digit-logical-shift-right
(%bignum-ref bignum i
)
1124 (%ashl
(%bignum-ref bignum
(1+ i
)) n-bits
))))))
1126 ;;; FIXNUM is assumed to be non-zero and the result of the shift should be a bignum
1127 (defun bignum-ashift-left-fixnum (fixnum count
)
1128 (declare (bignum-length count
)
1130 (multiple-value-bind (right-zero-digits remaining
)
1131 (truncate count digit-size
)
1132 (let* ((right-half (ldb (byte digit-size
0)
1133 (ash fixnum remaining
)))
1135 (logbitp (1- digit-size
) right-half
))
1136 (left-half (ash fixnum
1137 (- remaining digit-size
)))
1138 ;; Even if the left-half is 0 or -1 it might need to be sign
1139 ;; extended based on the left-most bit of the right-half
1140 (left-half-p (if sign-bit-p
1143 (length (+ right-zero-digits
1144 (if left-half-p
2 1)))
1145 (result (%allocate-bignum length
)))
1146 (setf (%bignum-ref result right-zero-digits
) right-half
)
1148 (setf (%bignum-ref result
(1+ right-zero-digits
))
1149 (ldb (byte digit-size
0) left-half
)))
1152 ;;;; relational operators
1154 ;;; This compares two bignums returning -1, 0, or 1, depending on
1155 ;;; whether a is less than, equal to, or greater than b.
1156 (declaim (ftype (function (bignum bignum
) (integer -
1 1)) bignum-compare
))
1157 (defun bignum-compare (a b
)
1158 (declare (type bignum a b
))
1159 (let* ((len-a (%bignum-length a
))
1160 (len-b (%bignum-length b
))
1161 (a-plusp (%bignum-0-or-plusp a len-a
))
1162 (b-plusp (%bignum-0-or-plusp b len-b
)))
1163 (declare (type bignum-length len-a len-b
))
1164 (cond ((not (eq a-plusp b-plusp
))
1167 (do ((i (1- len-a
) (1- i
)))
1169 (declare (type bignum-index i
))
1170 (let ((a-digit (%bignum-ref a i
))
1171 (b-digit (%bignum-ref b i
)))
1172 (declare (type bignum-element-type a-digit b-digit
))
1173 (when (> a-digit b-digit
)
1175 (when (> b-digit a-digit
)
1177 (when (zerop i
) (return 0))))
1180 (t (if a-plusp -
1 1)))))
1182 ;;;; float conversion
1184 ;;; Make a single or double float with the specified significand,
1185 ;;; exponent and sign.
1186 (defun single-float-from-bits (bits exp plusp
)
1187 (declare (fixnum exp
))
1188 (declare (optimize #-sb-xc-host
(inhibit-warnings 3)))
1190 sb
!vm
:single-float-exponent-byte
1191 (logandc2 (logand #xffffffff
1192 (%bignum-ref bits
1))
1193 sb
!vm
:single-float-hidden-bit
))))
1197 (logior res
(ash -
1 sb
!vm
:float-sign-shift
))))))
1198 (defun double-float-from-bits (bits exp plusp
)
1199 (declare (fixnum exp
))
1200 (declare (optimize #-sb-xc-host
(inhibit-warnings 3)))
1202 sb
!vm
:double-float-exponent-byte
1203 (logandc2 (ecase sb
!vm
::n-word-bits
1204 (32 (%bignum-ref bits
2))
1205 (64 (ash (%bignum-ref bits
1) -
32)))
1206 sb
!vm
:double-float-hidden-bit
)))
1207 (lo (logand #xffffffff
(%bignum-ref bits
1))))
1208 (make-double-float (if plusp
1210 (logior hi
(ash -
1 sb
!vm
:float-sign-shift
)))
1212 #!+(and long-float x86
)
1213 (defun long-float-from-bits (bits exp plusp
)
1214 (declare (fixnum exp
))
1215 (declare (optimize #-sb-xc-host
(inhibit-warnings 3)))
1219 (logior exp
(ash 1 15)))
1220 (%bignum-ref bits
2)
1221 (%bignum-ref bits
1)))
1223 ;;; Convert Bignum to a float in the specified Format, rounding to the best
1225 (defun bignum-to-float (bignum format
)
1226 (let* ((plusp (bignum-plus-p bignum
))
1227 (x (if plusp bignum
(negate-bignum bignum
)))
1228 (len (bignum-integer-length x
))
1229 (digits (float-format-digits format
))
1230 (keep (+ digits digit-size
))
1231 (shift (- keep len
))
1232 (shifted (if (minusp shift
)
1233 (bignum-ashift-right x
(- shift
))
1234 (bignum-ashift-left x shift
)))
1235 (low (%bignum-ref shifted
0))
1236 (round-bit (ash 1 (1- digit-size
))))
1237 (declare (type bignum-length len digits keep
) (fixnum shift
))
1238 (labels ((round-up ()
1239 (let ((rounded (add-bignums shifted round-bit
)))
1240 (if (> (integer-length rounded
) keep
)
1241 (float-from-bits (bignum-ashift-right rounded
1)
1243 (float-from-bits rounded len
))))
1244 (float-from-bits (bits len
)
1245 (declare (type bignum-length len
))
1248 (single-float-from-bits
1250 (check-exponent len sb
!vm
:single-float-bias
1251 sb
!vm
:single-float-normal-exponent-max
)
1254 (double-float-from-bits
1256 (check-exponent len sb
!vm
:double-float-bias
1257 sb
!vm
:double-float-normal-exponent-max
)
1261 (long-float-from-bits
1263 (check-exponent len sb
!vm
:long-float-bias
1264 sb
!vm
:long-float-normal-exponent-max
)
1266 (check-exponent (exp bias max
)
1267 (declare (type bignum-length len
))
1268 (let ((exp (+ exp bias
)))
1270 ;; Why a SIMPLE-TYPE-ERROR? Well, this is mainly
1271 ;; called by COERCE, which requires an error of
1272 ;; TYPE-ERROR if the conversion can't happen
1273 ;; (except in certain circumstances when we are
1274 ;; coercing to a FUNCTION) -- CSR, 2002-09-18
1275 (error 'simple-type-error
1276 :format-control
"Too large to be represented as a ~S:~% ~S"
1277 :format-arguments
(list format x
)
1278 :expected-type format
1283 ;; Round down if round bit is 0.
1284 ((not (logtest round-bit low
))
1285 (float-from-bits shifted len
))
1286 ;; If only round bit is set, then round to even.
1287 ((and (= low round-bit
)
1288 (dotimes (i (- (%bignum-length x
) (ceiling keep digit-size
))
1290 (unless (zerop (%bignum-ref x i
)) (return nil
))))
1291 (let ((next (%bignum-ref shifted
1)))
1294 (float-from-bits shifted len
))))
1295 ;; Otherwise, round up.
1299 ;;;; integer length and logbitp/logcount
1301 (defun bignum-buffer-integer-length (bignum len
)
1302 (declare (type bignum bignum
))
1303 (let* ((len-1 (1- len
))
1304 (digit (%bignum-ref bignum len-1
)))
1305 (declare (type bignum-length len len-1
)
1306 (type bignum-element-type digit
))
1307 (+ (integer-length (%fixnum-digit-with-correct-sign digit
))
1308 (* len-1 digit-size
))))
1310 (defun bignum-integer-length (bignum)
1311 (declare (type bignum bignum
))
1312 (bignum-buffer-integer-length bignum
(%bignum-length bignum
)))
1314 (defun bignum-logbitp (index bignum
)
1315 (declare (type bignum bignum
)
1316 (type bignum-index index
))
1317 (let ((len (%bignum-length bignum
)))
1318 (declare (type bignum-length len
))
1319 (multiple-value-bind (word-index bit-index
)
1320 (floor index digit-size
)
1321 (if (>= word-index len
)
1322 (not (bignum-plus-p bignum
))
1323 (logbitp bit-index
(%bignum-ref bignum word-index
))))))
1325 (defun bignum-logcount (bignum)
1326 (declare (type bignum bignum
)
1328 (declare (muffle-conditions compiler-note
)) ; returns lispobj, so what.
1329 (let ((length (%bignum-length bignum
))
1331 (declare (type bignum-length length
)
1333 (do ((index 0 (1+ index
)))
1335 (if (%bignum-0-or-plusp bignum length
)
1337 (- (* length digit-size
) result
)))
1338 (let ((digit (%bignum-ref bignum index
)))
1339 (declare (type bignum-element-type digit
))
1340 (incf result
(logcount digit
))))))
1342 ;;;; logical operations
1346 (defun bignum-logical-not (a)
1347 (declare (type bignum a
))
1348 (let* ((len (%bignum-length a
))
1349 (res (%allocate-bignum len
)))
1350 (declare (type bignum-length len
))
1351 (dotimes (i len res
)
1352 (declare (type bignum-index i
))
1353 (setf (%bignum-ref res i
) (%lognot
(%bignum-ref a i
))))))
1357 (defun bignum-logical-and (a b
)
1358 (declare (type bignum a b
))
1359 (let* ((len-a (%bignum-length a
))
1360 (len-b (%bignum-length b
))
1361 (a-plusp (%bignum-0-or-plusp a len-a
))
1362 (b-plusp (%bignum-0-or-plusp b len-b
)))
1363 (declare (type bignum-length len-a len-b
))
1367 (logand-shorter-positive a len-a b
(%allocate-bignum len-a
))
1368 (logand-shorter-negative a len-a b len-b
(%allocate-bignum len-b
))))
1371 (logand-shorter-positive b len-b a
(%allocate-bignum len-b
))
1372 (logand-shorter-negative b len-b a len-a
(%allocate-bignum len-a
))))
1373 (t (logand-shorter-positive a len-a b
(%allocate-bignum len-a
))))))
1375 ;;; This takes a shorter bignum, a and len-a, that is positive. Because this
1376 ;;; is AND, we don't care about any bits longer than a's since its infinite 0
1377 ;;; sign bits will mask the other bits out of b. The result is len-a big.
1378 (defun logand-shorter-positive (a len-a b res
)
1379 (declare (type bignum a b res
)
1380 (type bignum-length len-a
))
1382 (declare (type bignum-index i
))
1383 (setf (%bignum-ref res i
)
1384 (%logand
(%bignum-ref a i
) (%bignum-ref b i
))))
1385 (%normalize-bignum res len-a
))
1387 ;;; This takes a shorter bignum, a and len-a, that is negative. Because this
1388 ;;; is AND, we just copy any bits longer than a's since its infinite 1 sign
1389 ;;; bits will include any bits from b. The result is len-b big.
1390 (defun logand-shorter-negative (a len-a b len-b res
)
1391 (declare (type bignum a b res
)
1392 (type bignum-length len-a len-b
))
1394 (declare (type bignum-index i
))
1395 (setf (%bignum-ref res i
)
1396 (%logand
(%bignum-ref a i
) (%bignum-ref b i
))))
1397 (do ((i len-a
(1+ i
)))
1399 (declare (type bignum-index i
))
1400 (setf (%bignum-ref res i
) (%bignum-ref b i
)))
1401 (%normalize-bignum res len-b
))
1405 (defun bignum-logical-ior (a b
)
1406 (declare (type bignum a b
))
1407 (let* ((len-a (%bignum-length a
))
1408 (len-b (%bignum-length b
))
1409 (a-plusp (%bignum-0-or-plusp a len-a
))
1410 (b-plusp (%bignum-0-or-plusp b len-b
)))
1411 (declare (type bignum-length len-a len-b
))
1415 (logior-shorter-positive a len-a b len-b
(%allocate-bignum len-b
))
1416 (logior-shorter-negative a len-a b len-b
(%allocate-bignum len-b
))))
1419 (logior-shorter-positive b len-b a len-a
(%allocate-bignum len-a
))
1420 (logior-shorter-negative b len-b a len-a
(%allocate-bignum len-a
))))
1421 (t (logior-shorter-positive a len-a b len-b
(%allocate-bignum len-a
))))))
1423 ;;; This takes a shorter bignum, a and len-a, that is positive. Because this
1424 ;;; is IOR, we don't care about any bits longer than a's since its infinite
1425 ;;; 0 sign bits will mask the other bits out of b out to len-b. The result
1427 (defun logior-shorter-positive (a len-a b len-b res
)
1428 (declare (type bignum a b res
)
1429 (type bignum-length len-a len-b
))
1431 (declare (type bignum-index i
))
1432 (setf (%bignum-ref res i
)
1433 (%logior
(%bignum-ref a i
) (%bignum-ref b i
))))
1434 (do ((i len-a
(1+ i
)))
1436 (declare (type bignum-index i
))
1437 (setf (%bignum-ref res i
) (%bignum-ref b i
)))
1438 (%normalize-bignum res len-b
))
1440 ;;; This takes a shorter bignum, a and len-a, that is negative. Because this
1441 ;;; is IOR, we just copy any bits longer than a's since its infinite 1 sign
1442 ;;; bits will include any bits from b. The result is len-b long.
1443 (defun logior-shorter-negative (a len-a b len-b res
)
1444 (declare (type bignum a b res
)
1445 (type bignum-length len-a len-b
))
1447 (declare (type bignum-index i
))
1448 (setf (%bignum-ref res i
)
1449 (%logior
(%bignum-ref a i
) (%bignum-ref b i
))))
1450 (do ((i len-a
(1+ i
))
1451 (sign (%sign-digit a len-a
)))
1453 (declare (type bignum-index i
))
1454 (setf (%bignum-ref res i
) sign
))
1455 (%normalize-bignum res len-b
))
1459 (defun bignum-logical-xor (a b
)
1460 (declare (type bignum a b
))
1461 (let ((len-a (%bignum-length a
))
1462 (len-b (%bignum-length b
)))
1463 (declare (type bignum-length len-a len-b
))
1465 (bignum-logical-xor-aux a len-a b len-b
(%allocate-bignum len-b
))
1466 (bignum-logical-xor-aux b len-b a len-a
(%allocate-bignum len-a
)))))
1468 ;;; This takes the shorter of two bignums in a and len-a. Res is len-b
1469 ;;; long. Do the XOR.
1470 (defun bignum-logical-xor-aux (a len-a b len-b res
)
1471 (declare (type bignum a b res
)
1472 (type bignum-length len-a len-b
))
1474 (declare (type bignum-index i
))
1475 (setf (%bignum-ref res i
)
1476 (%logxor
(%bignum-ref a i
) (%bignum-ref b i
))))
1477 (do ((i len-a
(1+ i
))
1478 (sign (%sign-digit a len-a
)))
1480 (declare (type bignum-index i
))
1481 (setf (%bignum-ref res i
) (%logxor sign
(%bignum-ref b i
))))
1482 (%normalize-bignum res len-b
))
1484 ;;;; There used to be a bunch of code to implement "efficient" versions of LDB
1485 ;;;; and DPB here. But it apparently was never used, so it's been deleted.
1486 ;;;; --njf, 2007-02-04
1488 ;; This could be used by way of a transform, though for now it's specifically
1489 ;; a helper for %LDB in the limited case that it recognizes as non-consing.
1490 (defun ldb-bignum=>fixnum
(byte-size byte-pos bignum
)
1491 (declare (type (integer 0 #.sb
!vm
:n-positive-fixnum-bits
) byte-size
)
1492 (type bit-index byte-pos
))
1493 (multiple-value-bind (word-index bit-index
) (floor byte-pos digit-size
)
1494 (let ((n-digits (%bignum-length bignum
)))
1495 (cond ((>= word-index n-digits
) ; load from the infinitely extended sign word
1496 (ldb (byte byte-size
0) (%sign-digit bignum n-digits
)))
1497 ((<= (+ bit-index byte-size
) digit-size
) ; contained in one word
1498 ;; This case takes care of byte-size = 0 also.
1499 (ldb (byte byte-size bit-index
) (%bignum-ref bignum word-index
)))
1501 ;; At least one bit is obtained from each of two words,
1502 ;; and not more than two words.
1503 (let* ((low-part-size
1504 (truly-the (integer 1 #.
(1- sb
!vm
:n-positive-fixnum-bits
))
1505 (- digit-size bit-index
)))
1507 (truly-the (integer 1 #.
(1- sb
!vm
:n-positive-fixnum-bits
))
1508 (- byte-size low-part-size
))))
1509 (logior (truly-the (and fixnum unsigned-byte
) ; high part
1510 (let ((word-index (1+ word-index
)))
1511 (if (< word-index n-digits
) ; next word exists
1512 (ash (ldb (byte high-part-size
0)
1513 (%bignum-ref bignum word-index
))
1515 (mask-field (byte high-part-size low-part-size
)
1516 (%sign-digit bignum n-digits
)))))
1517 (ldb (byte low-part-size bit-index
) ; low part
1518 (%bignum-ref bignum word-index
)))))))))
1522 ;;; This is the original sketch of the algorithm from which I implemented this
1523 ;;; TRUNCATE, assuming both operands are bignums. I should modify this to work
1524 ;;; with the documentation on my functions, as a general introduction. I've
1525 ;;; left this here just in case someone needs it in the future. Don't look at
1526 ;;; this unless reading the functions' comments leaves you at a loss. Remember
1527 ;;; this comes from Knuth, so the book might give you the right general
1532 ;;; If X's magnitude is less than Y's, then result is 0 with remainder X.
1534 ;;; Make x and y positive, copying x if it is already positive.
1536 ;;; Shift y left until there's a 1 in the 30'th bit (most significant, non-sign
1538 ;;; Just do most sig digit to determine how much to shift whole number.
1539 ;;; Shift x this much too.
1540 ;;; Remember this initial shift count.
1542 ;;; Allocate q to be len-x minus len-y quantity plus 1.
1544 ;;; i = last digit of x.
1545 ;;; k = last digit of q.
1549 ;;; j = last digit of y.
1552 ;;; if x[i] = y[j] then g = (1- (ash 1 digit-size))
1553 ;;; else g = x[i]x[i-1]/y[j].
1556 ;;; %UNSIGNED-MULTIPLY returns b and c defined below.
1557 ;;; a = x[i-1] - (logand (* g y[j]) #xFFFFFFFF).
1558 ;;; Use %UNSIGNED-MULTIPLY taking low-order result.
1559 ;;; b = (logand (ash (* g y[j-1]) (- digit-size)) (1- (ash 1 digit-size))).
1560 ;;; c = (logand (* g y[j-1]) (1- (ash 1 digit-size))).
1562 ;;; if a > b, guess is too high
1563 ;;; g = g - 1; go back to "check guess".
1564 ;;; if a = b and c > x[i-2], guess is too high
1565 ;;; g = g - 1; go back to "check guess".
1566 ;;; GUESS IS 32-BIT NUMBER, SO USE THING TO KEEP IN SPECIAL REGISTER
1567 ;;; SAME FOR A, B, AND C.
1569 ;;; Subtract g * y from x[i - len-y+1]..x[i]. See paper for doing this in step.
1570 ;;; If x[i] < 0, guess is screwed up.
1571 ;;; negative g, then add 1
1572 ;;; zero or positive g, then subtract 1
1573 ;;; AND add y back into x[len-y+1..i].
1579 ;;; If k>=0, goto LOOP.
1581 ;;; Now quotient is good, but remainder is not.
1582 ;;; Shift x right by saved initial left shifting count.
1584 ;;; Check quotient and remainder signs.
1585 ;;; x pos y pos --> q pos r pos
1586 ;;; x pos y neg --> q neg r pos
1587 ;;; x neg y pos --> q neg r neg
1588 ;;; x neg y neg --> q pos r neg
1590 ;;; Normalize quotient and remainder. Cons result if necessary.
1593 ;;; This used to be split into multiple functions, which shared state
1594 ;;; in special variables *TRUNCATE-X* and *TRUNCATE-Y*. Having so many
1595 ;;; special variable accesses in tight inner loops was having a large
1596 ;;; effect on performance, so the helper functions have now been
1597 ;;; refactored into local functions and the special variables into
1598 ;;; lexicals. There was also a lot of boxing and unboxing of
1599 ;;; (UNSIGNED-BYTE 32)'s going on, which this refactoring
1600 ;;; eliminated. This improves the performance on some CL-BENCH tests
1601 ;;; by up to 50%, which is probably signigicant enough to justify the
1602 ;;; reduction in readability that was introduced. --JES, 2004-08-07
1603 (defun bignum-truncate (x y
)
1604 (declare (type bignum x y
))
1605 (declare (muffle-conditions compiler-note
)) ; returns lispobj, so what.
1606 (let (truncate-x truncate-y
)
1608 ;;; Divide X by Y when Y is a single bignum digit. BIGNUM-TRUNCATE
1609 ;;; fixes up the quotient and remainder with respect to sign and
1612 ;;; We don't have to worry about shifting Y to make its most
1613 ;;; significant digit sufficiently large for %BIGFLOOR to return
1614 ;;; digit-size quantities for the q-digit and r-digit. If Y is
1615 ;;; a single digit bignum, it is already large enough for
1616 ;;; %BIGFLOOR. That is, it has some bits on pretty high in the
1618 ((bignum-truncate-single-digit (x len-x y
)
1619 (declare (type bignum-length len-x
))
1620 (let ((y (%bignum-ref y
0)))
1621 (declare (type bignum-element-type y
))
1622 (if (not (logtest y
(1- y
)))
1623 ;; Y is a power of two.
1624 ;; SHIFT-RIGHT-UNALIGNED won't do the right thing
1625 ;; with a shift count of 0 or -1, so special case this.
1627 (error 'division-by-zero
:operation
'truncate
1628 :operands
(list x y
)))
1630 ;; We could probably get away with (VALUES X 0)
1631 ;; here, but it's not clear that some of the
1632 ;; normalization logic further down would avoid
1633 ;; mutilating X. Just go ahead and cons, consing's
1635 (values (copy-bignum x len-x
) 0))
1637 (let ((n-bits (1- (integer-length y
))))
1639 (shift-right-unaligned x
0 n-bits len-x
1641 (setf (%bignum-ref res j
)
1642 (%ashr
(%bignum-ref x i
) n-bits
))
1645 (logand (%bignum-ref x
0) (1- y
))))))
1646 (do ((i (1- len-x
) (1- i
))
1647 (q (%allocate-bignum len-x
))
1650 (let ((rem (%allocate-bignum
1)))
1651 (setf (%bignum-ref rem
0) r
)
1653 (declare (type bignum-element-type r
))
1654 (multiple-value-bind (q-digit r-digit
)
1655 (%bigfloor r
(%bignum-ref x i
) y
)
1656 (declare (type bignum-element-type q-digit r-digit
))
1657 (setf (%bignum-ref q i
) q-digit
)
1658 (setf r r-digit
))))))
1659 ;;; This returns a guess for the next division step. Y1 is the
1660 ;;; highest y digit, and y2 is the second to highest y
1661 ;;; digit. The x... variables are the three highest x digits
1662 ;;; for the next division step.
1664 ;;; From Knuth, our guess is either all ones or x-i and x-i-1
1665 ;;; divided by y1, depending on whether x-i and y1 are the
1666 ;;; same. We test this guess by determining whether guess*y2
1667 ;;; is greater than the three high digits of x minus guess*y1
1668 ;;; shifted left one digit:
1669 ;;; ------------------------------
1670 ;;; | x-i | x-i-1 | x-i-2 |
1671 ;;; ------------------------------
1672 ;;; ------------------------------
1673 ;;; - | g*y1 high | g*y1 low | 0 |
1674 ;;; ------------------------------
1675 ;;; ... < guess*y2 ???
1676 ;;; If guess*y2 is greater, then we decrement our guess by one
1677 ;;; and try again. This returns a guess that is either
1678 ;;; correct or one too large.
1679 (bignum-truncate-guess (y1 y2 x-i x-i-1 x-i-2
)
1680 (declare (type bignum-element-type y1 y2 x-i x-i-1 x-i-2
))
1681 (let ((guess (if (= x-i y1
)
1683 (%bigfloor x-i x-i-1 y1
))))
1684 (declare (type bignum-element-type guess
))
1686 (multiple-value-bind (high-guess*y1 low-guess
*y1
)
1687 (%multiply guess y1
)
1688 (declare (type bignum-element-type low-guess
*y1
1690 (multiple-value-bind (high-guess*y2 low-guess
*y2
)
1691 (%multiply guess y2
)
1692 (declare (type bignum-element-type high-guess
*y2
1694 (multiple-value-bind (middle-digit borrow
)
1695 (%subtract-with-borrow x-i-1 low-guess
*y1
1)
1696 (declare (type bignum-element-type middle-digit
)
1698 ;; Supplying borrow of 1 means there was no
1699 ;; borrow, and we know x-i-2 minus 0 requires
1701 (let ((high-digit (%subtract-with-borrow x-i
1704 (declare (type bignum-element-type high-digit
))
1705 (if (and (= high-digit
0)
1706 (or (> high-guess
*y2 middle-digit
)
1707 (and (= middle-digit high-guess
*y2
)
1708 (> low-guess
*y2 x-i-2
))))
1709 (setf guess
(%subtract-with-borrow guess
1 1))
1710 (return guess
)))))))))
1711 ;;; Divide TRUNCATE-X by TRUNCATE-Y, returning the quotient
1712 ;;; and destructively modifying TRUNCATE-X so that it holds
1715 ;;; LEN-X and LEN-Y tell us how much of the buffers we care about.
1717 ;;; TRUNCATE-X definitely has at least three digits, and it has one
1718 ;;; more than TRUNCATE-Y. This keeps i, i-1, i-2, and low-x-digit
1719 ;;; happy. Thanks to SHIFT-AND-STORE-TRUNCATE-BUFFERS.
1720 (return-quotient-leaving-remainder (len-x len-y
)
1721 (declare (type bignum-length len-x len-y
))
1722 (let* ((len-q (- len-x len-y
))
1723 ;; Add one for extra sign digit in case high bit is on.
1724 (q (%allocate-bignum
(1+ len-q
)))
1726 (y1 (%bignum-ref truncate-y
(1- len-y
)))
1727 (y2 (%bignum-ref truncate-y
(- len-y
2)))
1731 (low-x-digit (- i len-y
)))
1732 (declare (type bignum-length len-q
)
1733 (type bignum-index k i i-1 i-2 low-x-digit
)
1734 (type bignum-element-type y1 y2
))
1736 (setf (%bignum-ref q k
)
1737 (try-bignum-truncate-guess
1738 ;; This modifies TRUNCATE-X. Must access
1739 ;; elements each pass.
1740 (bignum-truncate-guess y1 y2
1741 (%bignum-ref truncate-x i
)
1742 (%bignum-ref truncate-x i-1
)
1743 (%bignum-ref truncate-x i-2
))
1745 (cond ((zerop k
) (return))
1748 (shiftf i i-1 i-2
(1- i-2
)))))
1750 ;;; This takes a digit guess, multiplies it by TRUNCATE-Y for a
1751 ;;; result one greater in length than LEN-Y, and subtracts this result
1752 ;;; from TRUNCATE-X. LOW-X-DIGIT is the first digit of X to start
1753 ;;; the subtraction, and we know X is long enough to subtract a LEN-Y
1754 ;;; plus one length bignum from it. Next we check the result of the
1755 ;;; subtraction, and if the high digit in X became negative, then our
1756 ;;; guess was one too big. In this case, return one less than GUESS
1757 ;;; passed in, and add one value of Y back into X to account for
1758 ;;; subtracting one too many. Knuth shows that the guess is wrong on
1759 ;;; the order of 3/b, where b is the base (2 to the digit-size power)
1760 ;;; -- pretty rarely.
1761 (try-bignum-truncate-guess (guess len-y low-x-digit
)
1762 (declare (type bignum-index low-x-digit
)
1763 (type bignum-length len-y
)
1764 (type bignum-element-type guess
))
1765 (let ((carry-digit 0)
1768 (declare (type bignum-element-type carry-digit
)
1769 (type bignum-index i
)
1771 ;; Multiply guess and divisor, subtracting from dividend
1774 (multiple-value-bind (high-digit low-digit
)
1775 (%multiply-and-add guess
1776 (%bignum-ref truncate-y j
)
1778 (declare (type bignum-element-type high-digit low-digit
))
1779 (setf carry-digit high-digit
)
1780 (multiple-value-bind (x temp-borrow
)
1781 (%subtract-with-borrow
(%bignum-ref truncate-x i
)
1784 (declare (type bignum-element-type x
)
1785 (fixnum temp-borrow
))
1786 (setf (%bignum-ref truncate-x i
) x
)
1787 (setf borrow temp-borrow
)))
1789 (setf (%bignum-ref truncate-x i
)
1790 (%subtract-with-borrow
(%bignum-ref truncate-x i
)
1791 carry-digit borrow
))
1792 ;; See whether guess is off by one, adding one
1793 ;; Y back in if necessary.
1794 (cond ((%digit-0-or-plusp
(%bignum-ref truncate-x i
))
1797 ;; If subtraction has negative result, add one
1798 ;; divisor value back in. The guess was one too
1799 ;; large in magnitude.
1800 (let ((i low-x-digit
)
1803 (multiple-value-bind (v k
)
1804 (%add-with-carry
(%bignum-ref truncate-y j
)
1805 (%bignum-ref truncate-x i
)
1807 (declare (type bignum-element-type v
))
1808 (setf (%bignum-ref truncate-x i
) v
)
1811 (setf (%bignum-ref truncate-x i
)
1812 (%add-with-carry
(%bignum-ref truncate-x i
)
1814 (%subtract-with-borrow guess
1 1)))))
1815 ;;; This returns the amount to shift y to place a one in the
1816 ;;; second highest bit. Y must be positive. If the last digit
1817 ;;; of y is zero, then y has a one in the previous digit's
1818 ;;; sign bit, so we know it will take one less than digit-size
1819 ;;; to get a one where we want. Otherwise, we count how many
1820 ;;; right shifts it takes to get zero; subtracting this value
1821 ;;; from digit-size tells us how many high zeros there are
1822 ;;; which is one more than the shift amount sought.
1824 ;;; Note: This is exactly the same as one less than the
1825 ;;; integer-length of the last digit subtracted from the
1828 ;;; We shift y to make it sufficiently large that doing the
1829 ;;; 2*digit-size by digit-size %BIGFLOOR calls ensures the quotient and
1830 ;;; remainder fit in digit-size.
1831 (shift-y-for-truncate (y)
1832 (let* ((len (%bignum-length y
))
1833 (last (%bignum-ref y
(1- len
))))
1834 (declare (type bignum-length len
)
1835 (type bignum-element-type last
))
1836 (- digit-size
(integer-length last
) 1)))
1837 ;;; Stores two bignums into the truncation bignum buffers,
1838 ;;; shifting them on the way in. This assumes x and y are
1839 ;;; positive and at least two in length, and it assumes
1840 ;;; truncate-x and truncate-y are one digit longer than x and
1842 (shift-and-store-truncate-buffers (x len-x y len-y shift
)
1843 (declare (type bignum-length len-x len-y
)
1844 (type (integer 0 (#.digit-size
)) shift
))
1845 (cond ((zerop shift
)
1846 (bignum-replace truncate-x x
:end1 len-x
)
1847 (bignum-replace truncate-y y
:end1 len-y
))
1849 (bignum-ashift-left-unaligned x
0 shift
(1+ len-x
)
1851 (bignum-ashift-left-unaligned y
0 shift
(1+ len-y
)
1852 truncate-y
))))) ;; LABELS
1853 ;;; Divide X by Y returning the quotient and remainder. In the
1854 ;;; general case, we shift Y to set up for the algorithm, and we
1855 ;;; use two buffers to save consing intermediate values. X gets
1856 ;;; destructively modified to become the remainder, and we have
1857 ;;; to shift it to account for the initial Y shift. After we
1858 ;;; multiple bind q and r, we first fix up the signs and then
1859 ;;; return the normalized results.
1860 (let* ((x-plusp (bignum-plus-p x
))
1861 (y-plusp (bignum-plus-p y
))
1862 (x (if x-plusp x
(negate-bignum x nil
)))
1863 (y (if y-plusp y
(negate-bignum y nil
)))
1864 (len-x (%bignum-length x
))
1865 (len-y (%bignum-length y
)))
1866 (multiple-value-bind (q r
)
1868 (bignum-truncate-single-digit x len-x y
))
1869 ((plusp (bignum-compare y x
))
1870 (let ((res (%allocate-bignum len-x
)))
1872 (setf (%bignum-ref res i
) (%bignum-ref x i
)))
1875 (let ((len-x+1 (1+ len-x
)))
1876 (setf truncate-x
(%allocate-bignum len-x
+1))
1877 (setf truncate-y
(%allocate-bignum
(1+ len-y
)))
1878 (let ((y-shift (shift-y-for-truncate y
)))
1879 (shift-and-store-truncate-buffers x len-x y
1881 (values (return-quotient-leaving-remainder len-x
+1
1883 ;; Now that RETURN-QUOTIENT-LEAVING-REMAINDER
1884 ;; has executed, we just tidy up the remainder
1885 ;; (in TRUNCATE-X) and return it.
1888 (let ((res (%allocate-bignum len-y
)))
1889 (declare (type bignum res
))
1890 (bignum-replace res truncate-x
:end2 len-y
)
1891 (%normalize-bignum res len-y
)))
1893 (shift-right-unaligned
1894 truncate-x
0 y-shift len-y
1896 (setf (%bignum-ref res j
)
1897 (%ashr
(%bignum-ref truncate-x i
)
1899 (%normalize-bignum res res-len
))
1901 (let ((quotient (cond ((eq x-plusp y-plusp
) q
)
1902 ((typep q
'fixnum
) (the fixnum
(- q
)))
1903 (t (negate-bignum-in-place q
))))
1904 (rem (cond (x-plusp r
)
1905 ((typep r
'fixnum
) (the fixnum
(- r
)))
1906 (t (negate-bignum-in-place r
)))))
1907 (values (if (typep quotient
'fixnum
)
1909 (%normalize-bignum quotient
(%bignum-length quotient
)))
1910 (if (typep rem
'fixnum
)
1912 (%normalize-bignum rem
(%bignum-length rem
))))))))))
1915 ;;;; There used to be a pile of code for implementing division for bignum digits
1916 ;;;; for machines that don't have a 2*digit-size by digit-size divide instruction.
1917 ;;;; This happens to be most machines, but all the SBCL ports seem to be content
1918 ;;;; to implement SB-BIGNUM:%BIGFLOOR as a VOP rather than using the code here.
1919 ;;;; So it's been deleted. --njf, 2007-02-04
1921 ;;;; general utilities
1923 ;;; Internal in-place operations use this to fixup remaining digits in the
1924 ;;; incoming data, such as in-place shifting. This is basically the same as
1925 ;;; the first form in %NORMALIZE-BIGNUM, but we return the length of the buffer
1926 ;;; instead of shrinking the bignum.
1927 #!-sb-fluid
(declaim (maybe-inline %normalize-bignum-buffer
))
1928 (defun %normalize-bignum-buffer
(result len
)
1929 (declare (type bignum result
)
1930 (type bignum-length len
))
1932 (do ((next-digit (%bignum-ref result
(- len
2))
1933 (%bignum-ref result
(- len
2)))
1934 (sign-digit (%bignum-ref result
(1- len
)) next-digit
))
1935 ((not (zerop (logxor sign-digit
(%ashr next-digit
(1- digit-size
))))))
1937 (setf (%bignum-ref result len
) 0)
1942 ;;; This drops the last digit if it is unnecessary sign information. It repeats
1943 ;;; this as needed, possibly ending with a fixnum. If the resulting length from
1944 ;;; shrinking is one, see whether our one word is a fixnum. Shift the possible
1945 ;;; fixnum bits completely out of the word, and compare this with shifting the
1946 ;;; sign bit all the way through. If the bits are all 1's or 0's in both words,
1947 ;;; then there are just sign bits between the fixnum bits and the sign bit. If
1948 ;;; we do have a fixnum, shift it over for the two low-tag bits.
1949 (defun %normalize-bignum
(result len
)
1950 (declare (type bignum result
)
1951 (type bignum-length len
)
1952 (muffle-conditions compiler-note
)
1953 #!-sb-fluid
(inline %normalize-bignum-buffer
))
1954 (let ((newlen (%normalize-bignum-buffer result len
)))
1955 (declare (type bignum-length newlen
))
1956 (unless (= newlen len
)
1957 (%bignum-set-length result newlen
))
1959 (let ((digit (%bignum-ref result
0)))
1960 (if (= (%ashr digit sb
!vm
:n-positive-fixnum-bits
)
1961 (%ashr digit
(1- digit-size
)))
1962 (%fixnum-digit-with-correct-sign digit
)
1966 ;;; This drops the last digit if it is unnecessary sign information. It
1967 ;;; repeats this as needed, possibly ending with a fixnum magnitude but never
1968 ;;; returning a fixnum.
1969 (defun %mostly-normalize-bignum
(result len
)
1970 (declare (type bignum result
)
1971 (type bignum-length len
)
1972 #!-sb-fluid
(inline %normalize-bignum-buffer
))
1973 (let ((newlen (%normalize-bignum-buffer result len
)))
1974 (declare (type bignum-length newlen
))
1975 (unless (= newlen len
)
1976 (%bignum-set-length result newlen
))
1981 ;;; the bignum case of the SXHASH function
1982 (defun sxhash-bignum (x)
1983 (let ((result 316495330))
1984 (declare (type fixnum result
))
1985 (dotimes (i (%bignum-length x
))
1986 (declare (type index i
))
1987 (let ((xi (%bignum-ref x i
)))
1989 (logand most-positive-fixnum