6129b4e403f8a9028018f1ace493ea572078cdf7
[sbcl.git] / src / code / bignum.lisp
1 ;;;; code to implement bignum support
2
3 ;;;; This software is part of the SBCL system. See the README file for
4 ;;;; more information.
5 ;;;;
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.
11
12 (in-package "SB!BIGNUM")
13 \f
14 ;;;; notes
15
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-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))
33
34 ;;; The following interfaces will either be assembler routines or code
35 ;;; sequences expanded into the code as basic bignum operations:
36 ;;;    General:
37 ;;;       %BIGNUM-LENGTH
38 ;;;       %ALLOCATE-BIGNUM
39 ;;;       %BIGNUM-REF
40 ;;;       %NORMALIZE-BIGNUM
41 ;;;       %BIGNUM-SET-LENGTH
42 ;;;       %FIXNUM-DIGIT-WITH-CORRECT-SIGN
43 ;;;       %SIGN-DIGIT
44 ;;;       %ASHR
45 ;;;       %ASHL
46 ;;;       %BIGNUM-0-OR-PLUSP
47 ;;;       %DIGIT-LOGICAL-SHIFT-RIGHT
48 ;;;    General (May not exist when done due to sole use in %-routines.)
49 ;;;       %DIGIT-0-OR-PLUSP
50 ;;;    Addition:
51 ;;;       %ADD-WITH-CARRY
52 ;;;    Subtraction:
53 ;;;       %SUBTRACT-WITH-BORROW
54 ;;;    Multiplication
55 ;;;       %MULTIPLY
56 ;;;    Negation
57 ;;;       %LOGNOT
58 ;;;    Shifting (in place)
59 ;;;       %NORMALIZE-BIGNUM-BUFFER
60 ;;;    GCD/Relational operators:
61 ;;;       %DIGIT-COMPARE
62 ;;;       %DIGIT-GREATER
63 ;;;    Relational operators:
64 ;;;       %LOGAND
65 ;;;       %LOGIOR
66 ;;;       %LOGXOR
67 ;;;    LDB
68 ;;;       %FIXNUM-TO-DIGIT
69 ;;;    TRUNCATE
70 ;;;       %FLOOR
71 ;;;
72 ;;; Note: The floating routines know about the float representation.
73 ;;;
74 ;;; PROBLEM 1:
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 32-bit
77 ;;; registers and number stack slots. I bind locals to these values, and I
78 ;;; use function on them -- ZEROP, ASH, etc.
79 ;;;
80 ;;; PROBLEM 2:
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.
85 ;;;
86 ;;; To do:
87 ;;;    fixnums
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.
96 ;;;       addition
97 ;;;       subtraction (twice)
98 ;;;       multiply
99 ;;;       GCD
100 ;;;    Write MASK-FIELD and DEPOSIT-FIELD in terms of logical operations.
101 ;;;    DIVIDE
102 ;;;       IF (/ x y) with bignums:
103 ;;;       do the truncate, and if rem is 0, return quotient.
104 ;;;       if rem is non-0
105 ;;;          gcd of x and y.
106 ;;;          "truncate" each by gcd, ignoring remainder 0.
107 ;;;          form ratio of each result, bottom is positive.
108 \f
109 ;;;; What's a bignum?
110
111 (defconstant digit-size 32)
112
113 (defconstant maximum-bignum-length (1- (ash 1 (- 32 sb!vm:n-widetag-bits))))
114 \f
115 ;;;; internal inline routines
116
117 ;;; %ALLOCATE-BIGNUM must zero all elements.
118 (defun %allocate-bignum (length)
119   (declare (type bignum-index length))
120   (%allocate-bignum length))
121
122 ;;; Extract the length of the bignum.
123 (defun %bignum-length (bignum)
124   (declare (type bignum-type bignum))
125   (%bignum-length bignum))
126
127 ;;; %BIGNUM-REF needs to access bignums as obviously as possible, and it needs
128 ;;; to be able to return 32 bits somewhere no one looks for real objects.
129 (defun %bignum-ref (bignum i)
130   (declare (type bignum-type bignum)
131            (type bignum-index i))
132   (%bignum-ref bignum i))
133 (defun %bignum-set (bignum i value)
134   (declare (type bignum-type bignum)
135            (type bignum-index i)
136            (type bignum-element-type value))
137   (%bignum-set bignum i value))
138
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)))
143
144 #!-sb-fluid (declaim (inline %bignum-0-or-plusp))
145 (defun %bignum-0-or-plusp (bignum len)
146   (declare (type bignum-type bignum)
147            (type bignum-index len))
148   (%digit-0-or-plusp (%bignum-ref bignum (1- len))))
149
150 ;;; This should be in assembler, and should not cons intermediate results. It
151 ;;; returns a 32bit digit and a carry resulting from adding together a, b, and
152 ;;; 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))
157
158 ;;; This should be in assembler, and should not cons intermediate results. It
159 ;;; returns a 32bit digit and a borrow resulting from subtracting b from a, and
160 ;;; subtracting a possible incoming borrow.
161 ;;;
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))
167
168 ;;; Multiply two digit-size (32-bit) numbers, returning a 64-bit result
169 ;;; split into two 32-bit quantities.
170 (defun %multiply (x y)
171   (declare (type bignum-element-type x y))
172   (%multiply x y))
173
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 at a
177 ;;; time possibly producing two carries) from adding these three digits get
178 ;;; added to the high digit from the multiply, producing the next carry digit.
179 ;;; Res-digit is optional since two uses of this primitive multiplies a single
180 ;;; digit bignum by a multiple digit bignum, and in this situation there is no
181 ;;; need for a result buffer accumulating partial results which is where the
182 ;;; res-digit comes from.
183 (defun %multiply-and-add (x-digit y-digit carry-in-digit
184                           &optional (res-digit 0))
185   (declare (type bignum-element-type x-digit y-digit res-digit carry-in-digit))
186   (%multiply-and-add x-digit y-digit carry-in-digit res-digit))
187
188 (defun %lognot (digit)
189   (declare (type bignum-element-type digit))
190   (%lognot digit))
191
192 ;;; Each of these does the 32-bit unsigned op.
193 #!-sb-fluid (declaim (inline %logand %logior %logxor))
194 (defun %logand (a b)
195   (declare (type bignum-element-type a b))
196   (logand a b))
197 (defun %logior (a b)
198   (declare (type bignum-element-type a b))
199   (logior a b))
200 (defun %logxor (a b)
201   (declare (type bignum-element-type a b))
202   (logxor a b))
203
204 ;;; This takes a fixnum and sets it up as an unsigned 32-bit quantity. In
205 ;;; the new system this will mean shifting it right two bits.
206 (defun %fixnum-to-digit (x)
207   (declare (fixnum x))
208   (logand x (1- (ash 1 digit-size))))
209
210 #!-32x16-divide
211 ;;; This takes three digits and returns the FLOOR'ed result of dividing the
212 ;;; first two as a 64-bit integer by the third.
213 ;;;
214 ;;; DO WEIRD let AND setq STUFF TO SLIME THE COMPILER INTO ALLOWING THE %FLOOR
215 ;;; TRANSFORM TO EXPAND INTO PSEUDO-ASSEMBLER FOR WHICH THE COMPILER CAN LATER
216 ;;; CORRECTLY ALLOCATE REGISTERS.
217 (defun %floor (a b c)
218   (let ((a a) (b b) (c c))
219     (declare (type bignum-element-type a b c))
220     (setq a a b b c c)
221     (%floor a b c)))
222
223 ;;; Convert the digit to a regular integer assuming that the digit is signed.
224 (defun %fixnum-digit-with-correct-sign (digit)
225   (declare (type bignum-element-type digit))
226   (if (logbitp (1- digit-size) digit)
227       (logior digit (ash -1 digit-size))
228       digit))
229
230 ;;; Do an arithmetic shift right of data even though bignum-element-type is
231 ;;; unsigned.
232 (defun %ashr (data count)
233   (declare (type bignum-element-type data)
234            (type (mod 32) count))
235   (%ashr data count))
236
237 ;;; This takes a 32-bit quantity and shifts it to the left, returning a 32-bit
238 ;;; quantity.
239 (defun %ashl (data count)
240   (declare (type bignum-element-type data)
241            (type (mod 32) count))
242   (%ashl data count))
243
244 ;;; Do an unsigned (logical) right shift of a digit by Count.
245 (defun %digit-logical-shift-right (data count)
246   (declare (type bignum-element-type data)
247            (type (mod 32) count))
248   (%digit-logical-shift-right data count))
249
250 ;;; Change the length of bignum to be newlen. Newlen must be the same or
251 ;;; smaller than the old length, and any elements beyond newlen must be zeroed.
252 (defun %bignum-set-length (bignum newlen)
253   (declare (type bignum-type bignum)
254            (type bignum-index newlen))
255   (%bignum-set-length bignum newlen))
256
257 ;;; This returns 0 or "-1" depending on whether the bignum is positive. This
258 ;;; is suitable for infinite sign extension to complete additions,
259 ;;; subtractions, negations, etc. This cannot return a -1 represented as
260 ;;; a negative fixnum since it would then have to low zeros.
261 #!-sb-fluid (declaim (inline %sign-digit))
262 (defun %sign-digit (bignum len)
263   (declare (type bignum-type bignum)
264            (type bignum-index len))
265   (%ashr (%bignum-ref bignum (1- len)) (1- digit-size)))
266
267 ;;; These take two 32 bit quantities and compare or contrast them without
268 ;;; wasting time with incorrect type checking.
269 #!-sb-fluid (declaim (inline %digit-compare %digit-greater))
270 (defun %digit-compare (x y)
271   (= x y))
272 (defun %digit-greater (x y)
273   (> x y))
274 \f
275 (declaim (optimize (speed 3) (safety 0)))
276 \f
277 ;;;; addition
278
279 (defun add-bignums (a b)
280   (declare (type bignum-type a b))
281   (let ((len-a (%bignum-length a))
282         (len-b (%bignum-length b)))
283     (declare (type bignum-index len-a len-b))
284     (multiple-value-bind (a len-a b len-b)
285         (if (> len-a len-b)
286             (values a len-a b len-b)
287             (values b len-b a len-a))
288       (declare (type bignum-type a b)
289                (type bignum-index len-a len-b))
290       (let* ((len-res (1+ len-a))
291              (res (%allocate-bignum len-res))
292              (carry 0))
293         (declare (type bignum-index len-res)
294                  (type bignum-type res)
295                  (type (mod 2) carry))
296         (dotimes (i len-b)
297           (declare (type bignum-index i))
298           (multiple-value-bind (v k)
299               (%add-with-carry (%bignum-ref a i) (%bignum-ref b i) carry)
300             (declare (type bignum-element-type v)
301                      (type (mod 2) k))
302             (setf (%bignum-ref res i) v)
303             (setf carry k)))
304         (if (/= len-a len-b)
305             (finish-add a res carry (%sign-digit b len-b) len-b len-a)
306             (setf (%bignum-ref res len-a)
307                   (%add-with-carry (%sign-digit a len-a)
308                                    (%sign-digit b len-b)
309                                    carry)))
310         (%normalize-bignum res len-res)))))
311
312 ;;; This takes the longer of two bignums and propagates the carry through its
313 ;;; remaining high order digits.
314 (defun finish-add (a res carry sign-digit-b start end)
315   (declare (type bignum-type a res)
316            (type (mod 2) carry)
317            (type bignum-element-type sign-digit-b)
318            (type bignum-index start end))
319   (do ((i start (1+ i)))
320       ((= i end)
321        (setf (%bignum-ref res end)
322              (%add-with-carry (%sign-digit a end) sign-digit-b carry)))
323     (declare (type bignum-index i))
324     (multiple-value-bind (v k)
325         (%add-with-carry (%bignum-ref a i) sign-digit-b carry)
326       (setf (%bignum-ref res i) v)
327       (setf carry k)))
328   (values))
329 \f
330 ;;;; subtraction
331
332 (eval-when (:compile-toplevel :execute)
333
334 ;;; This subtracts b from a plugging result into res. Return-fun is the
335 ;;; function to call that fixes up the result returning any useful values, such
336 ;;; as the result. This macro may evaluate its arguments more than once.
337 (sb!xc:defmacro subtract-bignum-loop (a len-a b len-b res len-res return-fun)
338   (let ((borrow (gensym))
339         (a-digit (gensym))
340         (a-sign (gensym))
341         (b-digit (gensym))
342         (b-sign (gensym))
343         (i (gensym))
344         (v (gensym))
345         (k (gensym)))
346     `(let* ((,borrow 1)
347             (,a-sign (%sign-digit ,a ,len-a))
348             (,b-sign (%sign-digit ,b ,len-b)))
349        (declare (type bignum-element-type ,a-sign ,b-sign))
350        (dotimes (,i ,len-res)
351          (declare (type bignum-index ,i))
352          (let ((,a-digit (if (< ,i ,len-a) (%bignum-ref ,a ,i) ,a-sign))
353                (,b-digit (if (< ,i ,len-b) (%bignum-ref ,b ,i) ,b-sign)))
354            (declare (type bignum-element-type ,a-digit ,b-digit))
355            (multiple-value-bind (,v ,k)
356                (%subtract-with-borrow ,a-digit ,b-digit ,borrow)
357              (setf (%bignum-ref ,res ,i) ,v)
358              (setf ,borrow ,k))))
359        (,return-fun ,res ,len-res))))
360
361 ) ;EVAL-WHEN
362
363 (defun subtract-bignum (a b)
364   (declare (type bignum-type a b))
365   (let* ((len-a (%bignum-length a))
366          (len-b (%bignum-length b))
367          (len-res (1+ (max len-a len-b)))
368          (res (%allocate-bignum len-res)))
369     (declare (type bignum-index len-a len-b len-res)) ;Test len-res for bounds?
370     (subtract-bignum-loop a len-a b len-b res len-res %normalize-bignum)))
371
372 ;;; Operations requiring a subtraction without the overhead of intermediate
373 ;;; results, such as GCD, use this. It assumes Result is big enough for the
374 ;;; result.
375 (defun subtract-bignum-buffers (a len-a b len-b result)
376   (declare (type bignum-type a b)
377            (type bignum-index len-a len-b))
378   (let ((len-res (max len-a len-b)))
379     (subtract-bignum-loop a len-a b len-b result len-res
380                           %normalize-bignum-buffer)))
381 \f
382 ;;;; multiplication
383
384 (defun multiply-bignums (a b)
385   (declare (type bignum-type a b))
386   (let* ((a-plusp (%bignum-0-or-plusp a (%bignum-length a)))
387          (b-plusp (%bignum-0-or-plusp b (%bignum-length b)))
388          (a (if a-plusp a (negate-bignum a)))
389          (b (if b-plusp b (negate-bignum b)))
390          (len-a (%bignum-length a))
391          (len-b (%bignum-length b))
392          (len-res (+ len-a len-b))
393          (res (%allocate-bignum len-res))
394          (negate-res (not (eq a-plusp b-plusp))))
395     (declare (type bignum-index len-a len-b len-res))
396     (dotimes (i len-a)
397       (declare (type bignum-index i))
398       (let ((carry-digit 0)
399             (x (%bignum-ref a i))
400             (k i))
401         (declare (type bignum-index k)
402                  (type bignum-element-type carry-digit x))
403         (dotimes (j len-b)
404           (multiple-value-bind (big-carry res-digit)
405               (%multiply-and-add x
406                                  (%bignum-ref b j)
407                                  (%bignum-ref res k)
408                                  carry-digit)
409             (declare (type bignum-element-type big-carry res-digit))
410             (setf (%bignum-ref res k) res-digit)
411             (setf carry-digit big-carry)
412             (incf k)))
413         (setf (%bignum-ref res k) carry-digit)))
414     (when negate-res (negate-bignum-in-place res))
415     (%normalize-bignum res len-res)))
416
417 (defun multiply-bignum-and-fixnum (bignum fixnum)
418   (declare (type bignum-type bignum) (type fixnum fixnum))
419   (let* ((bignum-plus-p (%bignum-0-or-plusp bignum (%bignum-length bignum)))
420          (fixnum-plus-p (not (minusp fixnum)))
421          (bignum (if bignum-plus-p bignum (negate-bignum bignum)))
422          (bignum-len (%bignum-length bignum))
423          (fixnum (if fixnum-plus-p fixnum (- fixnum)))
424          (result (%allocate-bignum (1+ bignum-len)))
425          (carry-digit 0))
426     (declare (type bignum-type bignum result)
427              (type bignum-index bignum-len)
428              (type bignum-element-type fixnum carry-digit))
429     (dotimes (index bignum-len)
430       (declare (type bignum-index index))
431       (multiple-value-bind (next-digit low)
432           (%multiply-and-add (%bignum-ref bignum index) fixnum carry-digit)
433         (declare (type bignum-element-type next-digit low))
434         (setf carry-digit next-digit)
435         (setf (%bignum-ref result index) low)))
436     (setf (%bignum-ref result bignum-len) carry-digit)
437     (unless (eq bignum-plus-p fixnum-plus-p)
438       (negate-bignum-in-place result))
439     (%normalize-bignum result (1+ bignum-len))))
440
441 (defun multiply-fixnums (a b)
442   (declare (fixnum a b))
443   (let* ((a-minusp (minusp a))
444          (b-minusp (minusp b)))
445     (multiple-value-bind (high low)
446         (%multiply (if a-minusp (- a) a)
447                    (if b-minusp (- b) b))
448       (declare (type bignum-element-type high low))
449       (if (and (zerop high)
450                (%digit-0-or-plusp low))
451           (let ((low (sb!ext:truly-the (unsigned-byte 31)
452                                        (%fixnum-digit-with-correct-sign low))))
453             (if (eq a-minusp b-minusp)
454                 low
455                 (- low)))
456           (let ((res (%allocate-bignum 2)))
457             (%bignum-set res 0 low)
458             (%bignum-set res 1 high)
459             (unless (eq a-minusp b-minusp) (negate-bignum-in-place res))
460             (%normalize-bignum res 2))))))
461 \f
462 ;;;; BIGNUM-REPLACE and WITH-BIGNUM-BUFFERS
463
464 (eval-when (:compile-toplevel :execute)
465
466 (sb!xc:defmacro bignum-replace (dest
467                                 src
468                                 &key
469                                 (start1 '0)
470                                 end1
471                                 (start2 '0)
472                                 end2
473                                 from-end)
474   (sb!int:once-only ((n-dest dest)
475                      (n-src src))
476     (let ((n-start1 (gensym))
477           (n-end1 (gensym))
478           (n-start2 (gensym))
479           (n-end2 (gensym))
480           (i1 (gensym))
481           (i2 (gensym))
482           (end1 (or end1 `(%bignum-length ,n-dest)))
483           (end2 (or end2 `(%bignum-length ,n-src))))
484       (if from-end
485           `(let ((,n-start1 ,start1)
486                  (,n-start2 ,start2))
487              (do ((,i1 (1- ,end1) (1- ,i1))
488                   (,i2 (1- ,end2) (1- ,i2)))
489                  ((or (< ,i1 ,n-start1) (< ,i2 ,n-start2)))
490                (declare (fixnum ,i1 ,i2))
491                (%bignum-set ,n-dest ,i1
492                             (%bignum-ref ,n-src ,i2))))
493           `(let ((,n-end1 ,end1)
494                  (,n-end2 ,end2))
495              (do ((,i1 ,start1 (1+ ,i1))
496                   (,i2 ,start2 (1+ ,i2)))
497                  ((or (>= ,i1 ,n-end1) (>= ,i2 ,n-end2)))
498                (declare (type bignum-index ,i1 ,i2))
499                (%bignum-set ,n-dest ,i1
500                             (%bignum-ref ,n-src ,i2))))))))
501
502 (sb!xc:defmacro with-bignum-buffers (specs &body body)
503   #!+sb-doc
504   "WITH-BIGNUM-BUFFERS ({(var size [init])}*) Form*"
505   (sb!int:collect ((binds)
506                    (inits))
507     (dolist (spec specs)
508       (let ((name (first spec))
509             (size (second spec)))
510         (binds `(,name (%allocate-bignum ,size)))
511         (let ((init (third spec)))
512           (when init
513             (inits `(bignum-replace ,name ,init))))))
514     `(let* ,(binds)
515        ,@(inits)
516        ,@body)))
517
518 ) ;EVAL-WHEN
519 \f
520 ;;;; GCD
521
522 (defun bignum-gcd (a b)
523   (declare (type bignum-type a b))
524   (let* ((a (if (%bignum-0-or-plusp a (%bignum-length a))
525                 a
526                 (negate-bignum a nil)))
527          (b (if (%bignum-0-or-plusp b (%bignum-length b))
528                 b
529                 (negate-bignum b nil)))
530          (len-a (%bignum-length a))
531          (len-b (%bignum-length b)))
532       (declare (type bignum-index len-a len-b))
533     (with-bignum-buffers ((a-buffer len-a a)
534                           (b-buffer len-b b)
535                           (res-buffer (max len-a len-b)))
536       (let* ((factors-of-two
537               (bignum-factors-of-two a-buffer len-a
538                                      b-buffer len-b))
539              (len-a (make-gcd-bignum-odd
540                      a-buffer
541                      (bignum-buffer-ashift-right a-buffer len-a
542                                                  factors-of-two)))
543              (len-b (make-gcd-bignum-odd
544                      b-buffer
545                      (bignum-buffer-ashift-right b-buffer len-b
546                                                  factors-of-two))))
547         (declare (type bignum-index len-a len-b))
548         (let ((x a-buffer)
549               (len-x len-a)
550               (y b-buffer)
551               (len-y len-b)
552               (z res-buffer))
553           (loop
554             (multiple-value-bind (u v len-v r len-r)
555                 (bignum-gcd-order-and-subtract x len-x y len-y z)
556               (declare (type bignum-index len-v len-r))
557               (when (and (= len-r 1) (zerop (%bignum-ref r 0)))
558                 (if (zerop factors-of-two)
559                     (let ((ret (%allocate-bignum len-v)))
560                       (dotimes (i len-v)
561                         (setf (%bignum-ref ret i) (%bignum-ref v i)))
562                       (return (%normalize-bignum ret len-v)))
563                     (return (bignum-ashift-left v factors-of-two len-v))))
564               (setf x v  len-x len-v)
565               (setf y r  len-y (make-gcd-bignum-odd r len-r))
566               (setf z u))))))))
567
568 (defun bignum-gcd-order-and-subtract (a len-a b len-b res)
569   (declare (type bignum-index len-a len-b) (type bignum-type a b))
570   (cond ((= len-a len-b)
571          (do ((i (1- len-a) (1- i)))
572              ((= i -1)
573               (setf (%bignum-ref res 0) 0)
574               (values a b len-b res 1))
575            (let ((a-digit (%bignum-ref a i))
576                  (b-digit (%bignum-ref b i)))
577              (cond ((%digit-compare a-digit b-digit))
578                    ((%digit-greater a-digit b-digit)
579                     (return
580                      (values a b len-b res
581                              (subtract-bignum-buffers a len-a b len-b res))))
582                    (t
583                     (return
584                      (values b a len-a res
585                              (subtract-bignum-buffers b len-b
586                                                       a len-a
587                                                       res))))))))
588         ((> len-a len-b)
589          (values a b len-b res
590                  (subtract-bignum-buffers a len-a b len-b res)))
591         (t
592          (values b a len-a res
593                  (subtract-bignum-buffers b len-b a len-a res)))))
594
595 (defun make-gcd-bignum-odd (a len-a)
596   (declare (type bignum-type a) (type bignum-index len-a))
597   (dotimes (index len-a)
598     (declare (type bignum-index index))
599     (do ((digit (%bignum-ref a index) (%ashr digit 1))
600          (increment 0 (1+ increment)))
601         ((zerop digit))
602       (declare (type (mod 32) increment))
603       (when (oddp digit)
604         (return-from make-gcd-bignum-odd
605                      (bignum-buffer-ashift-right a len-a
606                                                  (+ (* index digit-size)
607                                                     increment)))))))
608
609 (defun bignum-factors-of-two (a len-a b len-b)
610   (declare (type bignum-index len-a len-b) (type bignum-type a))
611   (do ((i 0 (1+ i))
612        (end (min len-a len-b)))
613       ((= i end) (error "Unexpected zero bignums?"))
614     (declare (type bignum-index i end))
615     (let ((or-digits (%logior (%bignum-ref a i) (%bignum-ref b i))))
616       (unless (zerop or-digits)
617         (return (do ((j 0 (1+ j))
618                      (or-digits or-digits (%ashr or-digits 1)))
619                     ((oddp or-digits) (+ (* i digit-size) j))
620                   (declare (type (mod 32) j))))))))
621 \f
622 ;;;; negation
623
624 (eval-when (:compile-toplevel :execute)
625
626 ;;; This negates bignum-len digits of bignum, storing the resulting digits into
627 ;;; result (possibly EQ to bignum) and returning whatever end-carry there is.
628 (sb!xc:defmacro bignum-negate-loop (bignum
629                                     bignum-len
630                                     &optional (result nil resultp))
631   (let ((carry (gensym))
632         (end (gensym))
633         (value (gensym))
634         (last (gensym)))
635     `(let* (,@(if (not resultp) `(,last))
636             (,carry
637              (multiple-value-bind (,value ,carry)
638                  (%add-with-carry (%lognot (%bignum-ref ,bignum 0)) 1 0)
639                ,(if resultp
640                     `(setf (%bignum-ref ,result 0) ,value)
641                     `(setf ,last ,value))
642                ,carry))
643             (i 1)
644             (,end ,bignum-len))
645        (declare (type bit ,carry)
646                 (type bignum-index i ,end))
647        (loop
648          (when (= i ,end) (return))
649          (multiple-value-bind (,value temp)
650              (%add-with-carry (%lognot (%bignum-ref ,bignum i)) 0 ,carry)
651            ,(if resultp
652                 `(setf (%bignum-ref ,result i) ,value)
653                 `(setf ,last ,value))
654            (setf ,carry temp))
655          (incf i))
656        ,(if resultp carry `(values ,carry ,last)))))
657
658 ) ; EVAL-WHEN
659
660 ;;; Fully-normalize is an internal optional. It cause this to always return
661 ;;; a bignum, without any extraneous digits, and it never returns a fixnum.
662 (defun negate-bignum (x &optional (fully-normalize t))
663   (declare (type bignum-type x))
664   (let* ((len-x (%bignum-length x))
665          (len-res (1+ len-x))
666          (res (%allocate-bignum len-res)))
667     (declare (type bignum-index len-x len-res)) ;Test len-res for range?
668     (let ((carry (bignum-negate-loop x len-x res)))
669       (setf (%bignum-ref res len-x)
670             (%add-with-carry (%lognot (%sign-digit x len-x)) 0 carry)))
671     (if fully-normalize
672         (%normalize-bignum res len-res)
673         (%mostly-normalize-bignum res len-res))))
674
675 ;;; This assumes bignum is positive; that is, the result of negating it will
676 ;;; stay in the provided allocated bignum.
677 (defun negate-bignum-in-place (bignum)
678   (bignum-negate-loop bignum (%bignum-length bignum) bignum)
679   bignum)
680 \f
681 ;;;; shifting
682
683 (defconstant all-ones-digit #xFFFFFFFF)
684
685 (eval-when (:compile-toplevel :execute)
686
687 ;;; This macro is used by BIGNUM-ASHIFT-RIGHT, BIGNUM-BUFFER-ASHIFT-RIGHT, and
688 ;;; BIGNUM-LDB-BIGNUM-RES. They supply a termination form that references
689 ;;; locals established by this form. Source is the source bignum. Start-digit
690 ;;; is the first digit in source from which we pull bits. Start-pos is the
691 ;;; first bit we want. Res-len-form is the form that computes the length of
692 ;;; the resulting bignum. Termination is a DO termination form with a test and
693 ;;; body. When result is supplied, it is the variable to which this binds a
694 ;;; newly allocated bignum.
695 ;;;
696 ;;; Given start-pos, 1-31 inclusively, of shift, we form the j'th resulting
697 ;;; digit from high bits of the i'th source digit and the start-pos number of
698 ;;; bits from the i+1'th source digit.
699 (sb!xc:defmacro shift-right-unaligned (source
700                                        start-digit
701                                        start-pos
702                                        res-len-form
703                                        termination
704                                        &optional result)
705   `(let* ((high-bits-in-first-digit (- digit-size ,start-pos))
706           (res-len ,res-len-form)
707           (res-len-1 (1- res-len))
708           ,@(if result `((,result (%allocate-bignum res-len)))))
709      (declare (type bignum-index res-len res-len-1))
710      (do ((i ,start-digit i+1)
711           (i+1 (1+ ,start-digit) (1+ i+1))
712           (j 0 (1+ j)))
713          ,termination
714        (declare (type bignum-index i i+1 j))
715        (setf (%bignum-ref ,(if result result source) j)
716              (%logior (%digit-logical-shift-right (%bignum-ref ,source i)
717                                                   ,start-pos)
718                       (%ashl (%bignum-ref ,source i+1)
719                              high-bits-in-first-digit))))))
720
721 ) ; EVAL-WHEN
722
723 ;;; First compute the number of whole digits to shift, shifting them by
724 ;;; skipping them when we start to pick up bits, and the number of bits to
725 ;;; shift the remaining digits into place. If the number of digits is greater
726 ;;; than the length of the bignum, then the result is either 0 or -1. If we
727 ;;; shift on a digit boundary (that is, n-bits is zero), then we just copy
728 ;;; digits. The last branch handles the general case which uses a macro that a
729 ;;; couple other routines use. The fifth argument to the macro references
730 ;;; locals established by the macro.
731 (defun bignum-ashift-right (bignum count)
732   (declare (type bignum-type bignum)
733            (type unsigned-byte count))
734   (let ((bignum-len (%bignum-length bignum)))
735     (declare (type bignum-index bignum-len))
736     (cond ((fixnump count)
737            (multiple-value-bind (digits n-bits) (truncate count digit-size)
738              (declare (type bignum-index digits))
739              (cond
740               ((>= digits bignum-len)
741                (if (%bignum-0-or-plusp bignum bignum-len) 0 -1))
742               ((zerop n-bits)
743                (bignum-ashift-right-digits bignum digits))
744               (t
745                (shift-right-unaligned bignum digits n-bits (- bignum-len digits)
746                                       ((= j res-len-1)
747                                        (setf (%bignum-ref res j)
748                                              (%ashr (%bignum-ref bignum i) n-bits))
749                                        (%normalize-bignum res res-len))
750                                       res)))))
751           ((> count bignum-len)
752            0)
753            ;; Since a FIXNUM should be big enough to address anything in
754            ;; memory, including arrays of bits, and since arrays of bits
755            ;; take up about the same space as corresponding fixnums, there
756            ;; should be no way that we fall through to this case: any shift
757            ;; right by a bignum should give zero. But let's check anyway:
758           (t (error "bignum overflow: can't shift right by ~S")))))
759
760 (defun bignum-ashift-right-digits (bignum digits)
761   (declare (type bignum-type bignum)
762            (type bignum-index digits))
763   (let* ((res-len (- (%bignum-length bignum) digits))
764          (res (%allocate-bignum res-len)))
765     (declare (type bignum-index res-len)
766              (type bignum-type res))
767     (bignum-replace res bignum :start2 digits)
768     (%normalize-bignum res res-len)))
769
770 ;;; GCD uses this for an in-place shifting operation. This is different enough
771 ;;; from BIGNUM-ASHIFT-RIGHT that it isn't worth folding the bodies into a
772 ;;; macro, but they share the basic algorithm. This routine foregoes a first
773 ;;; test for digits being greater than or equal to bignum-len since that will
774 ;;; never happen for its uses in GCD. We did fold the last branch into a macro
775 ;;; since it was duplicated a few times, and the fifth argument to it
776 ;;; references locals established by the macro.
777 (defun bignum-buffer-ashift-right (bignum bignum-len x)
778   (declare (type bignum-index bignum-len) (fixnum x))
779   (multiple-value-bind (digits n-bits) (truncate x digit-size)
780     (declare (type bignum-index digits))
781     (cond
782      ((zerop n-bits)
783       (let ((new-end (- bignum-len digits)))
784         (bignum-replace bignum bignum :end1 new-end :start2 digits
785                         :end2 bignum-len)
786         (%normalize-bignum-buffer bignum new-end)))
787      (t
788       (shift-right-unaligned bignum digits n-bits (- bignum-len digits)
789                              ((= j res-len-1)
790                               (setf (%bignum-ref bignum j)
791                                     (%ashr (%bignum-ref bignum i) n-bits))
792                               (%normalize-bignum-buffer bignum res-len)))))))
793
794 ;;; This handles shifting a bignum buffer to provide fresh bignum data for some
795 ;;; internal routines. We know bignum is safe when called with bignum-len.
796 ;;; First we compute the number of whole digits to shift, shifting them
797 ;;; starting to store farther along the result bignum. If we shift on a digit
798 ;;; boundary (that is, n-bits is zero), then we just copy digits. The last
799 ;;; branch handles the general case.
800 (defun bignum-ashift-left (bignum x &optional bignum-len)
801   (declare (type bignum-type bignum)
802            (type unsigned-byte x)
803            (type (or null bignum-index) bignum-len))
804   (if (fixnump x)
805     (multiple-value-bind (digits n-bits) (truncate x digit-size)
806       (let* ((bignum-len (or bignum-len (%bignum-length bignum)))
807              (res-len (+ digits bignum-len 1)))
808         (when (> res-len maximum-bignum-length)
809           (error "can't represent result of left shift"))
810         (if (zerop n-bits)
811           (bignum-ashift-left-digits bignum bignum-len digits)
812           (bignum-ashift-left-unaligned bignum digits n-bits res-len))))
813     ;; Left shift by a number too big to be represented as a fixnum
814     ;; would exceed our memory capacity, since a fixnum is big enough
815     ;; index any array, including a bit array.
816     (error "can't represent result of left shift")))
817
818 (defun bignum-ashift-left-digits (bignum bignum-len digits)
819   (declare (type bignum-index bignum-len digits))
820   (let* ((res-len (+ bignum-len digits))
821          (res (%allocate-bignum res-len)))
822     (declare (type bignum-index res-len))
823     (bignum-replace res bignum :start1 digits :end1 res-len :end2 bignum-len
824                     :from-end t)
825     res))
826
827 ;;; BIGNUM-TRUNCATE uses this to store into a bignum buffer by supplying res.
828 ;;; When res comes in non-nil, then this foregoes allocating a result, and it
829 ;;; normalizes the buffer instead of the would-be allocated result.
830 ;;;
831 ;;; We start storing into one digit higher than digits, storing a whole result
832 ;;; digit from parts of two contiguous digits from bignum. When the loop
833 ;;; finishes, we store the remaining bits from bignum's first digit in the
834 ;;; first non-zero result digit, digits. We also grab some left over high
835 ;;; bits from the last digit of bignum.
836 (defun bignum-ashift-left-unaligned (bignum digits n-bits res-len
837                                      &optional (res nil resp))
838   (declare (type bignum-index digits res-len)
839            (type (mod #.digit-size) n-bits))
840   (let* ((remaining-bits (- digit-size n-bits))
841          (res-len-1 (1- res-len))
842          (res (or res (%allocate-bignum res-len))))
843     (declare (type bignum-index res-len res-len-1))
844     (do ((i 0 i+1)
845          (i+1 1 (1+ i+1))
846          (j (1+ digits) (1+ j)))
847         ((= j res-len-1)
848          (setf (%bignum-ref res digits)
849                (%ashl (%bignum-ref bignum 0) n-bits))
850          (setf (%bignum-ref res j)
851                (%ashr (%bignum-ref bignum i) remaining-bits))
852          (if resp
853              (%normalize-bignum-buffer res res-len)
854              (%normalize-bignum res res-len)))
855       (declare (type bignum-index i i+1 j))
856       (setf (%bignum-ref res j)
857             (%logior (%digit-logical-shift-right (%bignum-ref bignum i)
858                                                  remaining-bits)
859                      (%ashl (%bignum-ref bignum i+1) n-bits))))))
860 \f
861 ;;;; relational operators
862
863 ;;; Return T iff bignum is positive.
864 (defun bignum-plus-p (bignum)
865   (declare (type bignum-type bignum))
866   (%bignum-0-or-plusp bignum (%bignum-length bignum)))
867
868 ;;; This compares two bignums returning -1, 0, or 1, depending on
869 ;;; whether a is less than, equal to, or greater than b.
870 (declaim (ftype (function (bignum bignum) (integer -1 1)) bignum-compare))
871 (defun bignum-compare (a b)
872   (declare (type bignum-type a b))
873   (let* ((len-a (%bignum-length a))
874          (len-b (%bignum-length b))
875          (a-plusp (%bignum-0-or-plusp a len-a))
876          (b-plusp (%bignum-0-or-plusp b len-b)))
877     (declare (type bignum-index len-a len-b))
878     (cond ((not (eq a-plusp b-plusp))
879            (if a-plusp 1 -1))
880           ((= len-a len-b)
881            (do ((i (1- len-a) (1- i)))
882                (())
883              (declare (type bignum-index i))
884              (let ((a-digit (%bignum-ref a i))
885                    (b-digit (%bignum-ref b i)))
886                (declare (type bignum-element-type a-digit b-digit))
887                (when (%digit-greater a-digit b-digit)
888                  (return 1))
889                (when (%digit-greater b-digit a-digit)
890                  (return -1)))
891              (when (zerop i) (return 0))))
892           ((> len-a len-b)
893            (if a-plusp 1 -1))
894           (t (if a-plusp -1 1)))))
895 \f
896 ;;;; float conversion
897
898 ;;; Make a single or double float with the specified significand,
899 ;;; exponent and sign.
900 (defun single-float-from-bits (bits exp plusp)
901   (declare (fixnum exp))
902   (declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
903   (let ((res (dpb exp
904                   sb!vm:single-float-exponent-byte
905                   (logandc2 (sb!ext:truly-the (unsigned-byte 31)
906                                               (%bignum-ref bits 1))
907                             sb!vm:single-float-hidden-bit))))
908     (make-single-float
909      (if plusp
910          res
911          (logior res (ash -1 sb!vm:float-sign-shift))))))
912 (defun double-float-from-bits (bits exp plusp)
913   (declare (fixnum exp))
914   (declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
915   (let ((hi (dpb exp
916                  sb!vm:double-float-exponent-byte
917                  (logandc2 (sb!ext:truly-the (unsigned-byte 31)
918                                              (%bignum-ref bits 2))
919                            sb!vm:double-float-hidden-bit))))
920     (make-double-float
921      (if plusp
922          hi
923          (logior hi (ash -1 sb!vm:float-sign-shift)))
924      (%bignum-ref bits 1))))
925 #!+(and long-float x86)
926 (defun long-float-from-bits (bits exp plusp)
927   (declare (fixnum exp))
928   (declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
929   (make-long-float
930    (if plusp
931        exp
932        (logior exp (ash 1 15)))
933    (%bignum-ref bits 2)
934    (%bignum-ref bits 1)))
935
936 ;;; Convert Bignum to a float in the specified Format, rounding to the best
937 ;;; approximation.
938 (defun bignum-to-float (bignum format)
939   (let* ((plusp (bignum-plus-p bignum))
940          (x (if plusp bignum (negate-bignum bignum)))
941          (len (bignum-integer-length x))
942          (digits (float-format-digits format))
943          (keep (+ digits digit-size))
944          (shift (- keep len))
945          (shifted (if (minusp shift)
946                       (bignum-ashift-right x (- shift))
947                       (bignum-ashift-left x shift)))
948          (low (%bignum-ref shifted 0))
949          (round-bit (ash 1 (1- digit-size))))
950     (declare (type bignum-index len digits keep) (fixnum shift))
951     (labels ((round-up ()
952                (let ((rounded (add-bignums shifted round-bit)))
953                  (if (> (integer-length rounded) keep)
954                      (float-from-bits (bignum-ashift-right rounded 1)
955                                       (1+ len))
956                      (float-from-bits rounded len))))
957              (float-from-bits (bits len)
958                (declare (type bignum-index len))
959                (ecase format
960                  (single-float
961                   (single-float-from-bits
962                    bits
963                    (check-exponent len sb!vm:single-float-bias
964                                    sb!vm:single-float-normal-exponent-max)
965                    plusp))
966                  (double-float
967                   (double-float-from-bits
968                    bits
969                    (check-exponent len sb!vm:double-float-bias
970                                    sb!vm:double-float-normal-exponent-max)
971                    plusp))
972                  #!+long-float
973                  (long-float
974                   (long-float-from-bits
975                    bits
976                    (check-exponent len sb!vm:long-float-bias
977                                    sb!vm:long-float-normal-exponent-max)
978                    plusp))))
979              (check-exponent (exp bias max)
980                (declare (type bignum-index len))
981                (let ((exp (+ exp bias)))
982                  (when (> exp max)
983                    (error "Too large to be represented as a ~S:~%  ~S"
984                           format x))
985                  exp)))
986
987     (cond
988      ;; Round down if round bit is 0.
989      ((zerop (logand round-bit low))
990       (float-from-bits shifted len))
991      ;; If only round bit is set, then round to even.
992      ((and (= low round-bit)
993            (dotimes (i (- (%bignum-length x) (ceiling keep digit-size))
994                        t)
995              (unless (zerop (%bignum-ref x i)) (return nil))))
996       (let ((next (%bignum-ref shifted 1)))
997         (if (oddp next)
998             (round-up)
999             (float-from-bits shifted len))))
1000      ;; Otherwise, round up.
1001      (t
1002       (round-up))))))
1003 \f
1004 ;;;; integer length and logcount
1005
1006 (defun bignum-integer-length (bignum)
1007   (declare (type bignum-type bignum))
1008   (let* ((len (%bignum-length bignum))
1009          (len-1 (1- len))
1010          (digit (%bignum-ref bignum len-1)))
1011     (declare (type bignum-index len len-1)
1012              (type bignum-element-type digit))
1013     (+ (integer-length (%fixnum-digit-with-correct-sign digit))
1014        (* len-1 digit-size))))
1015
1016 (defun bignum-logcount (bignum)
1017   (declare (type bignum-type bignum))
1018   (let* ((length (%bignum-length bignum))
1019          (plusp (%bignum-0-or-plusp bignum length))
1020          (result 0))
1021     (declare (type bignum-index length)
1022              (fixnum result))
1023     (do ((index 0 (1+ index)))
1024         ((= index length) result)
1025       (let ((digit (%bignum-ref bignum index)))
1026         (declare (type bignum-element-type digit))
1027         (incf result (logcount (if plusp digit (%lognot digit))))))))
1028 \f
1029 ;;;; logical operations
1030
1031 ;;;; NOT
1032
1033 (defun bignum-logical-not (a)
1034   (declare (type bignum-type a))
1035   (let* ((len (%bignum-length a))
1036          (res (%allocate-bignum len)))
1037     (declare (type bignum-index len))
1038     (dotimes (i len res)
1039       (declare (type bignum-index i))
1040       (setf (%bignum-ref res i) (%lognot (%bignum-ref a i))))))
1041
1042 ;;;; AND
1043
1044 (defun bignum-logical-and (a b)
1045   (declare (type bignum-type a b))
1046   (let* ((len-a (%bignum-length a))
1047          (len-b (%bignum-length b))
1048          (a-plusp (%bignum-0-or-plusp a len-a))
1049          (b-plusp (%bignum-0-or-plusp b len-b)))
1050     (declare (type bignum-index len-a len-b))
1051     (cond
1052      ((< len-a len-b)
1053       (if a-plusp
1054           (logand-shorter-positive a len-a b (%allocate-bignum len-a))
1055           (logand-shorter-negative a len-a b len-b (%allocate-bignum len-b))))
1056      ((< len-b len-a)
1057       (if b-plusp
1058           (logand-shorter-positive b len-b a (%allocate-bignum len-b))
1059           (logand-shorter-negative b len-b a len-a (%allocate-bignum len-a))))
1060      (t (logand-shorter-positive a len-a b (%allocate-bignum len-a))))))
1061
1062 ;;; This takes a shorter bignum, a and len-a, that is positive. Because this
1063 ;;; is AND, we don't care about any bits longer than a's since its infinite 0
1064 ;;; sign bits will mask the other bits out of b. The result is len-a big.
1065 (defun logand-shorter-positive (a len-a b res)
1066   (declare (type bignum-type a b res)
1067            (type bignum-index len-a))
1068   (dotimes (i len-a)
1069     (declare (type bignum-index i))
1070     (setf (%bignum-ref res i)
1071           (%logand (%bignum-ref a i) (%bignum-ref b i))))
1072   (%normalize-bignum res len-a))
1073
1074 ;;; This takes a shorter bignum, a and len-a, that is negative. Because this
1075 ;;; is AND, we just copy any bits longer than a's since its infinite 1 sign
1076 ;;; bits will include any bits from b. The result is len-b big.
1077 (defun logand-shorter-negative (a len-a b len-b res)
1078   (declare (type bignum-type a b res)
1079            (type bignum-index len-a len-b))
1080   (dotimes (i len-a)
1081     (declare (type bignum-index i))
1082     (setf (%bignum-ref res i)
1083           (%logand (%bignum-ref a i) (%bignum-ref b i))))
1084   (do ((i len-a (1+ i)))
1085       ((= i len-b))
1086     (declare (type bignum-index i))
1087     (setf (%bignum-ref res i) (%bignum-ref b i)))
1088   (%normalize-bignum res len-b))
1089
1090 ;;;; IOR
1091
1092 (defun bignum-logical-ior (a b)
1093   (declare (type bignum-type a b))
1094   (let* ((len-a (%bignum-length a))
1095          (len-b (%bignum-length b))
1096          (a-plusp (%bignum-0-or-plusp a len-a))
1097          (b-plusp (%bignum-0-or-plusp b len-b)))
1098     (declare (type bignum-index len-a len-b))
1099     (cond
1100      ((< len-a len-b)
1101       (if a-plusp
1102           (logior-shorter-positive a len-a b len-b (%allocate-bignum len-b))
1103           (logior-shorter-negative a len-a b len-b (%allocate-bignum len-b))))
1104      ((< len-b len-a)
1105       (if b-plusp
1106           (logior-shorter-positive b len-b a len-a (%allocate-bignum len-a))
1107           (logior-shorter-negative b len-b a len-a (%allocate-bignum len-a))))
1108      (t (logior-shorter-positive a len-a b len-b (%allocate-bignum len-a))))))
1109
1110 ;;; This takes a shorter bignum, a and len-a, that is positive. Because this
1111 ;;; is IOR, we don't care about any bits longer than a's since its infinite
1112 ;;; 0 sign bits will mask the other bits out of b out to len-b. The result
1113 ;;; is len-b long.
1114 (defun logior-shorter-positive (a len-a b len-b res)
1115   (declare (type bignum-type a b res)
1116            (type bignum-index len-a len-b))
1117   (dotimes (i len-a)
1118     (declare (type bignum-index i))
1119     (setf (%bignum-ref res i)
1120           (%logior (%bignum-ref a i) (%bignum-ref b i))))
1121   (do ((i len-a (1+ i)))
1122       ((= i len-b))
1123     (declare (type bignum-index i))
1124     (setf (%bignum-ref res i) (%bignum-ref b i)))
1125   (%normalize-bignum res len-b))
1126
1127 ;;; This takes a shorter bignum, a and len-a, that is negative. Because this
1128 ;;; is IOR, we just copy any bits longer than a's since its infinite 1 sign
1129 ;;; bits will include any bits from b. The result is len-b long.
1130 (defun logior-shorter-negative (a len-a b len-b res)
1131   (declare (type bignum-type a b res)
1132            (type bignum-index len-a len-b))
1133   (dotimes (i len-a)
1134     (declare (type bignum-index i))
1135     (setf (%bignum-ref res i)
1136           (%logior (%bignum-ref a i) (%bignum-ref b i))))
1137   (do ((i len-a (1+ i))
1138        (sign (%sign-digit a len-a)))
1139       ((= i len-b))
1140     (declare (type bignum-index i))
1141     (setf (%bignum-ref res i) sign))
1142   (%normalize-bignum res len-b))
1143
1144 ;;;; XOR
1145
1146 (defun bignum-logical-xor (a b)
1147   (declare (type bignum-type a b))
1148   (let ((len-a (%bignum-length a))
1149         (len-b (%bignum-length b)))
1150     (declare (type bignum-index len-a len-b))
1151     (if (< len-a len-b)
1152         (bignum-logical-xor-aux a len-a b len-b (%allocate-bignum len-b))
1153         (bignum-logical-xor-aux b len-b a len-a (%allocate-bignum len-a)))))
1154
1155 ;;; This takes the shorter of two bignums in a and len-a. Res is len-b
1156 ;;; long. Do the XOR.
1157 (defun bignum-logical-xor-aux (a len-a b len-b res)
1158   (declare (type bignum-type a b res)
1159            (type bignum-index len-a len-b))
1160   (dotimes (i len-a)
1161     (declare (type bignum-index i))
1162     (setf (%bignum-ref res i)
1163           (%logxor (%bignum-ref a i) (%bignum-ref b i))))
1164   (do ((i len-a (1+ i))
1165        (sign (%sign-digit a len-a)))
1166       ((= i len-b))
1167     (declare (type bignum-index i))
1168     (setf (%bignum-ref res i) (%logxor sign (%bignum-ref b i))))
1169   (%normalize-bignum res len-b))
1170 \f
1171 ;;;; LDB (load byte)
1172
1173 #|
1174 FOR NOW WE DON'T USE LDB OR DPB. WE USE SHIFTS AND MASKS IN NUMBERS.LISP WHICH
1175 IS LESS EFFICIENT BUT EASIER TO MAINTAIN. BILL SAYS THIS CODE CERTAINLY WORKS!
1176
1177 (defconstant maximum-fixnum-bits #!+ibm-rt-pc 27 #!-ibm-rt-pc 30)
1178
1179 (defun bignum-load-byte (byte bignum)
1180   (declare (type bignum-type bignum))
1181   (let ((byte-len (byte-size byte))
1182         (byte-pos (byte-position byte)))
1183     (if (< byte-len maximum-fixnum-bits)
1184         (bignum-ldb-fixnum-res bignum byte-len byte-pos)
1185         (bignum-ldb-bignum-res bignum byte-len byte-pos))))
1186
1187 ;;; This returns a fixnum result of loading a byte from a bignum. In order, we
1188 ;;; check for the following conditions:
1189 ;;;    Insufficient bignum digits to start loading a byte --
1190 ;;;       Return 0 or byte-len 1's depending on sign of bignum.
1191 ;;;    One bignum digit containing the whole byte spec --
1192 ;;;       Grab 'em, shift 'em, and mask out what we don't want.
1193 ;;;    Insufficient bignum digits to cover crossing a digit boundary --
1194 ;;;       Grab the available bits in the last digit, and or in whatever
1195 ;;;       virtual sign bits we need to return a full byte spec.
1196 ;;;    Else (we cross a digit boundary with all bits available) --
1197 ;;;       Make a couple masks, grab what we want, shift it around, and
1198 ;;;       LOGIOR it all together.
1199 ;;; Because (< maximum-fixnum-bits digit-size) and
1200 ;;;      (< byte-len maximum-fixnum-bits),
1201 ;;; we only cross one digit boundary if any.
1202 (defun bignum-ldb-fixnum-res (bignum byte-len byte-pos)
1203   (multiple-value-bind (skipped-digits pos) (truncate byte-pos digit-size)
1204     (let ((bignum-len (%bignum-length bignum))
1205           (s-digits+1 (1+ skipped-digits)))
1206       (declare (type bignum-index bignum-len s-digits+1))
1207       (if (>= skipped-digits bignum-len)
1208           (if (%bignum-0-or-plusp bignum bignum-len)
1209               0
1210               (%make-ones byte-len))
1211           (let ((end (+ pos byte-len)))
1212             (cond ((<= end digit-size)
1213                    (logand (ash (%bignum-ref bignum skipped-digits) (- pos))
1214                            ;; Must LOGAND after shift here.
1215                            (%make-ones byte-len)))
1216                   ((>= s-digits+1 bignum-len)
1217                    (let* ((available-bits (- digit-size pos))
1218                           (res (logand (ash (%bignum-ref bignum skipped-digits)
1219                                             (- pos))
1220                                        ;; LOGAND should be unnecessary here
1221                                        ;; with a logical right shift or a
1222                                        ;; correct unsigned-byte-32 one.
1223                                        (%make-ones available-bits))))
1224                      (if (%bignum-0-or-plusp bignum bignum-len)
1225                          res
1226                          (logior (%ashl (%make-ones (- end digit-size))
1227                                         available-bits)
1228                                  res))))
1229                   (t
1230                    (let* ((high-bits-in-first-digit (- digit-size pos))
1231                           (high-mask (%make-ones high-bits-in-first-digit))
1232                           (low-bits-in-next-digit (- end digit-size))
1233                           (low-mask (%make-ones low-bits-in-next-digit)))
1234                      (declare (type bignum-element-type high-mask low-mask))
1235                      (logior (%ashl (logand (%bignum-ref bignum s-digits+1)
1236                                             low-mask)
1237                                     high-bits-in-first-digit)
1238                              (logand (ash (%bignum-ref bignum skipped-digits)
1239                                           (- pos))
1240                                      ;; LOGAND should be unnecessary here with
1241                                      ;; a logical right shift or a correct
1242                                      ;; unsigned-byte-32 one.
1243                                      high-mask))))))))))
1244
1245 ;;; This returns a bignum result of loading a byte from a bignum. In order, we
1246 ;;; check for the following conditions:
1247 ;;;    Insufficient bignum digits to start loading a byte --
1248 ;;;    Byte-pos starting on a digit boundary --
1249 ;;;    Byte spec contained in one bignum digit --
1250 ;;;       Grab the bits we want and stick them in a single digit result.
1251 ;;;       Since we know byte-pos is non-zero here, we know our single digit
1252 ;;;       will have a zero high sign bit.
1253 ;;;    Else (unaligned multiple digits) --
1254 ;;;       This is like doing a shift right combined with either masking
1255 ;;;       out unwanted high bits from bignum or filling in virtual sign
1256 ;;;       bits if bignum had insufficient bits. We use SHIFT-RIGHT-ALIGNED
1257 ;;;       and reference lots of local variables this macro establishes.
1258 (defun bignum-ldb-bignum-res (bignum byte-len byte-pos)
1259   (multiple-value-bind (skipped-digits pos) (truncate byte-pos digit-size)
1260     (let ((bignum-len (%bignum-length bignum)))
1261       (declare (type bignum-index bignum-len))
1262       (cond
1263        ((>= skipped-digits bignum-len)
1264         (make-bignum-virtual-ldb-bits bignum bignum-len byte-len))
1265        ((zerop pos)
1266         (make-aligned-ldb-bignum bignum bignum-len byte-len skipped-digits))
1267        ((< (+ pos byte-len) digit-size)
1268         (let ((res (%allocate-bignum 1)))
1269           (setf (%bignum-ref res 0)
1270                 (logand (%ashr (%bignum-ref bignum skipped-digits) pos)
1271                         (%make-ones byte-len)))
1272           res))
1273        (t
1274         (make-unaligned-ldb-bignum bignum bignum-len
1275                                    byte-len skipped-digits pos))))))
1276
1277 ;;; This returns bits from bignum that don't physically exist. These are
1278 ;;; all zero or one depending on the sign of the bignum.
1279 (defun make-bignum-virtual-ldb-bits (bignum bignum-len byte-len)
1280   (if (%bignum-0-or-plusp bignum bignum-len)
1281       0
1282       (multiple-value-bind (res-len-1 extra) (truncate byte-len digit-size)
1283         (declare (type bignum-index res-len-1))
1284         (let* ((res-len (1+ res-len-1))
1285                (res (%allocate-bignum res-len)))
1286           (declare (type bignum-index res-len))
1287           (do ((j 0 (1+ j)))
1288               ((= j res-len-1)
1289                (setf (%bignum-ref res j) (%make-ones extra))
1290                (%normalize-bignum res res-len))
1291             (declare (type bignum-index j))
1292             (setf (%bignum-ref res j) all-ones-digit))))))
1293
1294 ;;; Since we are picking up aligned digits, we just copy the whole digits
1295 ;;; we want and fill in extra bits. We might have a byte-len that extends
1296 ;;; off the end of the bignum, so we may have to fill in extra 1's if the
1297 ;;; bignum is negative.
1298 (defun make-aligned-ldb-bignum (bignum bignum-len byte-len skipped-digits)
1299   (multiple-value-bind (res-len-1 extra) (truncate byte-len digit-size)
1300     (declare (type bignum-index res-len-1))
1301     (let* ((res-len (1+ res-len-1))
1302            (res (%allocate-bignum res-len)))
1303       (declare (type bignum-index res-len))
1304       (do ((i skipped-digits (1+ i))
1305            (j 0 (1+ j)))
1306           ((or (= j res-len-1) (= i bignum-len))
1307            (cond ((< i bignum-len)
1308                   (setf (%bignum-ref res j)
1309                         (logand (%bignum-ref bignum i)
1310                                 (the bignum-element-type (%make-ones extra)))))
1311                  ((%bignum-0-or-plusp bignum bignum-len))
1312                  (t
1313                   (do ((j j (1+ j)))
1314                       ((= j res-len-1)
1315                        (setf (%bignum-ref res j) (%make-ones extra)))
1316                     (setf (%bignum-ref res j) all-ones-digit))))
1317            (%normalize-bignum res res-len))
1318       (declare (type bignum-index i j))
1319       (setf (%bignum-ref res j) (%bignum-ref bignum i))))))
1320
1321 ;;; This grabs unaligned bignum bits from bignum assuming byte-len causes at
1322 ;;; least one digit boundary crossing. We use SHIFT-RIGHT-UNALIGNED referencing
1323 ;;; lots of local variables established by it.
1324 (defun make-unaligned-ldb-bignum (bignum
1325                                   bignum-len
1326                                   byte-len
1327                                   skipped-digits
1328                                   pos)
1329   (multiple-value-bind (res-len-1 extra) (truncate byte-len digit-size)
1330     (shift-right-unaligned
1331      bignum skipped-digits pos (1+ res-len-1)
1332      ((or (= j res-len-1) (= i+1 bignum-len))
1333       (cond ((= j res-len-1)
1334              (cond
1335               ((< extra high-bits-in-first-digit)
1336                (setf (%bignum-ref res j)
1337                      (logand (ash (%bignum-ref bignum i) minus-start-pos)
1338                              ;; Must LOGAND after shift here.
1339                              (%make-ones extra))))
1340               (t
1341                (setf (%bignum-ref res j)
1342                      (logand (ash (%bignum-ref bignum i) minus-start-pos)
1343                              ;; LOGAND should be unnecessary here with a logical
1344                              ;; right shift or a correct unsigned-byte-32 one.
1345                              high-mask))
1346                (when (%bignum-0-or-plusp bignum bignum-len)
1347                  (setf (%bignum-ref res j)
1348                        (logior (%bignum-ref res j)
1349                                (%ashl (%make-ones
1350                                        (- extra high-bits-in-first-digit))
1351                                       high-bits-in-first-digit)))))))
1352             (t
1353              (setf (%bignum-ref res j)
1354                    (logand (ash (%bignum-ref bignum i) minus-start-pos)
1355                            ;; LOGAND should be unnecessary here with a logical
1356                            ;; right shift or a correct unsigned-byte-32 one.
1357                            high-mask))
1358              (unless (%bignum-0-or-plusp bignum bignum-len)
1359                ;; Fill in upper half of this result digit with 1's.
1360                (setf (%bignum-ref res j)
1361                      (logior (%bignum-ref res j)
1362                              (%ashl low-mask high-bits-in-first-digit)))
1363                ;; Fill in any extra 1's we need to be byte-len long.
1364                (do ((j (1+ j) (1+ j)))
1365                    ((>= j res-len-1)
1366                     (setf (%bignum-ref res j) (%make-ones extra)))
1367                  (setf (%bignum-ref res j) all-ones-digit)))))
1368       (%normalize-bignum res res-len))
1369      res)))
1370 \f
1371 ;;;; DPB (deposit byte)
1372
1373 (defun bignum-deposit-byte (new-byte byte-spec bignum)
1374   (declare (type bignum-type bignum))
1375   (let* ((byte-len (byte-size byte-spec))
1376          (byte-pos (byte-position byte-spec))
1377          (bignum-len (%bignum-length bignum))
1378          (bignum-plusp (%bignum-0-or-plusp bignum bignum-len))
1379          (byte-end (+ byte-pos byte-len))
1380          (res-len (1+ (max (ceiling byte-end digit-size) bignum-len)))
1381          (res (%allocate-bignum res-len)))
1382     (declare (type bignum-index bignum-len res-len))
1383     ;; Fill in an extra sign digit in case we set what would otherwise be the
1384     ;; last digit's last bit. Normalize at the end in case this was
1385     ;; unnecessary.
1386     (unless bignum-plusp
1387       (setf (%bignum-ref res (1- res-len)) all-ones-digit))
1388     (multiple-value-bind (end-digit end-bits) (truncate byte-end digit-size)
1389       (declare (type bignum-index end-digit))
1390       ;; Fill in bits from bignum up to byte-pos.
1391       (multiple-value-bind (pos-digit pos-bits) (truncate byte-pos digit-size)
1392         (declare (type bignum-index pos-digit))
1393         (do ((i 0 (1+ i))
1394              (end (min pos-digit bignum-len)))
1395             ((= i end)
1396              (cond ((< i bignum-len)
1397                     (unless (zerop pos-bits)
1398                       (setf (%bignum-ref res i)
1399                             (logand (%bignum-ref bignum i)
1400                                     (%make-ones pos-bits)))))
1401                    (bignum-plusp)
1402                    (t
1403                     (do ((i i (1+ i)))
1404                         ((= i pos-digit)
1405                          (unless (zerop pos-bits)
1406                            (setf (%bignum-ref res i) (%make-ones pos-bits))))
1407                       (setf (%bignum-ref res i) all-ones-digit)))))
1408           (setf (%bignum-ref res i) (%bignum-ref bignum i)))
1409         ;; Fill in bits from new-byte.
1410         (if (typep new-byte 'fixnum)
1411             (deposit-fixnum-bits new-byte byte-len pos-digit pos-bits
1412                                  end-digit end-bits res)
1413             (deposit-bignum-bits new-byte byte-len pos-digit pos-bits
1414                                  end-digit end-bits res)))
1415       ;; Fill in remaining bits from bignum after byte-spec.
1416       (when (< end-digit bignum-len)
1417         (setf (%bignum-ref res end-digit)
1418               (logior (logand (%bignum-ref bignum end-digit)
1419                               (%ashl (%make-ones (- digit-size end-bits))
1420                                      end-bits))
1421                       ;; DEPOSIT-FIXNUM-BITS and DEPOSIT-BIGNUM-BITS only store
1422                       ;; bits from new-byte into res's end-digit element, so
1423                       ;; we don't need to mask out unwanted high bits.
1424                       (%bignum-ref res end-digit)))
1425         (do ((i (1+ end-digit) (1+ i)))
1426             ((= i bignum-len))
1427           (setf (%bignum-ref res i) (%bignum-ref bignum i)))))
1428     (%normalize-bignum res res-len)))
1429
1430 ;;; This starts at result's pos-digit skipping pos-bits, and it stores bits
1431 ;;; from new-byte, a fixnum, into result. It effectively stores byte-len
1432 ;;; number of bits, but never stores past end-digit and end-bits in result.
1433 ;;; The first branch fires when all the bits we want from new-byte are present;
1434 ;;; if byte-len crosses from the current result digit into the next, the last
1435 ;;; argument to DEPOSIT-FIXNUM-DIGIT is a mask for those bits. The second
1436 ;;; branch handles the need to grab more bits than the fixnum new-byte has, but
1437 ;;; new-byte is positive; therefore, any virtual bits are zero. The mask for
1438 ;;; bits that don't fit in the current result digit is simply the remaining
1439 ;;; bits in the bignum digit containing new-byte; we don't care if we store
1440 ;;; some extra in the next result digit since they will be zeros. The last
1441 ;;; branch handles the need to grab more bits than the fixnum new-byte has, but
1442 ;;; new-byte is negative; therefore, any virtual bits must be explicitly filled
1443 ;;; in as ones. We call DEPOSIT-FIXNUM-DIGIT to grab what bits actually exist
1444 ;;; and to fill in the current result digit.
1445 (defun deposit-fixnum-bits (new-byte byte-len pos-digit pos-bits
1446                             end-digit end-bits result)
1447   (declare (type bignum-index pos-digit end-digit))
1448   (let ((other-bits (- digit-size pos-bits))
1449         (new-byte-digit (%fixnum-to-digit new-byte)))
1450     (declare (type bignum-element-type new-byte-digit))
1451     (cond ((< byte-len maximum-fixnum-bits)
1452            (deposit-fixnum-digit new-byte-digit byte-len pos-digit pos-bits
1453                                  other-bits result
1454                                  (- byte-len other-bits)))
1455           ((or (plusp new-byte) (zerop new-byte))
1456            (deposit-fixnum-digit new-byte-digit byte-len pos-digit pos-bits
1457                                  other-bits result pos-bits))
1458           (t
1459            (multiple-value-bind (digit bits)
1460                (deposit-fixnum-digit new-byte-digit byte-len pos-digit pos-bits
1461                                      other-bits result
1462                                      (if (< (- byte-len other-bits) digit-size)
1463                                          (- byte-len other-bits)
1464                                          digit-size))
1465              (declare (type bignum-index digit))
1466              (cond ((< digit end-digit)
1467                     (setf (%bignum-ref result digit)
1468                           (logior (%bignum-ref result digit)
1469                                   (%ashl (%make-ones (- digit-size bits)) bits)))
1470                     (do ((i (1+ digit) (1+ i)))
1471                         ((= i end-digit)
1472                          (setf (%bignum-ref result i) (%make-ones end-bits)))
1473                       (setf (%bignum-ref result i) all-ones-digit)))
1474                    ((> digit end-digit))
1475                    ((< bits end-bits)
1476                     (setf (%bignum-ref result digit)
1477                           (logior (%bignum-ref result digit)
1478                                   (%ashl (%make-ones (- end-bits bits))
1479                                          bits))))))))))
1480
1481 ;;; This fills in the current result digit from new-byte-digit. The first case
1482 ;;; handles everything we want fitting in the current digit, and other-bits is
1483 ;;; the number of bits remaining to be filled in result's current digit. This
1484 ;;; number is digit-size minus pos-bits. The second branch handles filling in
1485 ;;; result's current digit, and it shoves the unused bits of new-byte-digit
1486 ;;; into the next result digit. This is correct regardless of new-byte-digit's
1487 ;;; sign. It returns the new current result digit and how many bits already
1488 ;;; filled in the result digit.
1489 (defun deposit-fixnum-digit (new-byte-digit byte-len pos-digit pos-bits
1490                              other-bits result next-digit-bits-needed)
1491   (declare (type bignum-index pos-digit)
1492            (type bignum-element-type new-byte-digit next-digit-mask))
1493   (cond ((<= byte-len other-bits)
1494          ;; Bits from new-byte fit in the current result digit.
1495          (setf (%bignum-ref result pos-digit)
1496                (logior (%bignum-ref result pos-digit)
1497                        (%ashl (logand new-byte-digit (%make-ones byte-len))
1498                               pos-bits)))
1499          (if (= byte-len other-bits)
1500              (values (1+ pos-digit) 0)
1501              (values pos-digit (+ byte-len pos-bits))))
1502         (t
1503          ;; Some of new-byte's bits go in current result digit.
1504          (setf (%bignum-ref result pos-digit)
1505                (logior (%bignum-ref result pos-digit)
1506                        (%ashl (logand new-byte-digit (%make-ones other-bits))
1507                               pos-bits)))
1508          (let ((pos-digit+1 (1+ pos-digit)))
1509            ;; The rest of new-byte's bits go in the next result digit.
1510            (setf (%bignum-ref result pos-digit+1)
1511                  (logand (ash new-byte-digit (- other-bits))
1512                          ;; Must LOGAND after shift here.
1513                          (%make-ones next-digit-bits-needed)))
1514            (if (= next-digit-bits-needed digit-size)
1515                (values (1+ pos-digit+1) 0)
1516                (values pos-digit+1 next-digit-bits-needed))))))
1517
1518 ;;; This starts at result's pos-digit skipping pos-bits, and it stores bits
1519 ;;; from new-byte, a bignum, into result. It effectively stores byte-len
1520 ;;; number of bits, but never stores past end-digit and end-bits in result.
1521 ;;; When handling a starting bit unaligned with a digit boundary, we check
1522 ;;; in the second branch for the byte spec fitting into the pos-digit element
1523 ;;; after after pos-bits; DEPOSIT-UNALIGNED-BIGNUM-BITS expects at least one
1524 ;;; digit boundary crossing.
1525 (defun deposit-bignum-bits (bignum-byte byte-len pos-digit pos-bits
1526                             end-digit end-bits result)
1527   (declare (type bignum-index pos-digit end-digit))
1528   (cond ((zerop pos-bits)
1529          (deposit-aligned-bignum-bits bignum-byte pos-digit end-digit end-bits
1530                                       result))
1531         ((or (= end-digit pos-digit)
1532              (and (= end-digit (1+ pos-digit))
1533                   (zerop end-bits)))
1534          (setf (%bignum-ref result pos-digit)
1535                (logior (%bignum-ref result pos-digit)
1536                        (%ashl (logand (%bignum-ref bignum-byte 0)
1537                                       (%make-ones byte-len))
1538                               pos-bits))))
1539         (t (deposit-unaligned-bignum-bits bignum-byte pos-digit pos-bits
1540                                           end-digit end-bits result))))
1541
1542 ;;; This deposits bits from bignum-byte into result starting at pos-digit and
1543 ;;; the zero'th bit. It effectively only stores bits to end-bits in the
1544 ;;; end-digit element of result. The loop termination code takes care of
1545 ;;; picking up the last digit's bits or filling in virtual negative sign bits.
1546 (defun deposit-aligned-bignum-bits (bignum-byte pos-digit end-digit end-bits
1547                                     result)
1548   (declare (type bignum-index pos-digit end-digit))
1549   (let* ((bignum-len (%bignum-length bignum-byte))
1550          (bignum-plusp (%bignum-0-or-plusp bignum-byte bignum-len)))
1551     (declare (type bignum-index bignum-len))
1552     (do ((i 0 (1+ i ))
1553          (j pos-digit (1+ j)))
1554         ((or (= j end-digit) (= i bignum-len))
1555          (cond ((= j end-digit)
1556                 (cond ((< i bignum-len)
1557                        (setf (%bignum-ref result j)
1558                              (logand (%bignum-ref bignum-byte i)
1559                                      (%make-ones end-bits))))
1560                       (bignum-plusp)
1561                       (t
1562                        (setf (%bignum-ref result j) (%make-ones end-bits)))))
1563                (bignum-plusp)
1564                (t
1565                 (do ((j j (1+ j)))
1566                     ((= j end-digit)
1567                      (setf (%bignum-ref result j) (%make-ones end-bits)))
1568                   (setf (%bignum-ref result j) all-ones-digit)))))
1569       (setf (%bignum-ref result j) (%bignum-ref bignum-byte i)))))
1570
1571 ;;; This assumes at least one digit crossing.
1572 (defun deposit-unaligned-bignum-bits (bignum-byte pos-digit pos-bits
1573                                       end-digit end-bits result)
1574   (declare (type bignum-index pos-digit end-digit))
1575   (let* ((bignum-len (%bignum-length bignum-byte))
1576          (bignum-plusp (%bignum-0-or-plusp bignum-byte bignum-len))
1577          (low-mask (%make-ones pos-bits))
1578          (bits-past-pos-bits (- digit-size pos-bits))
1579          (high-mask (%make-ones bits-past-pos-bits))
1580          (minus-high-bits (- bits-past-pos-bits)))
1581     (declare (type bignum-element-type low-mask high-mask)
1582              (type bignum-index bignum-len))
1583     (do ((i 0 (1+ i))
1584          (j pos-digit j+1)
1585          (j+1 (1+ pos-digit) (1+ j+1)))
1586         ((or (= j end-digit) (= i bignum-len))
1587          (cond
1588           ((= j end-digit)
1589            (setf (%bignum-ref result j)
1590                  (cond
1591                   ((>= pos-bits end-bits)
1592                    (logand (%bignum-ref result j) (%make-ones end-bits)))
1593                   ((< i bignum-len)
1594                    (logior (%bignum-ref result j)
1595                            (%ashl (logand (%bignum-ref bignum-byte i)
1596                                           (%make-ones (- end-bits pos-bits)))
1597                                   pos-bits)))
1598                   (bignum-plusp
1599                    (logand (%bignum-ref result j)
1600                            ;; 0's between pos-bits and end-bits positions.
1601                            (logior (%ashl (%make-ones (- digit-size end-bits))
1602                                           end-bits)
1603                                    low-mask)))
1604                   (t (logior (%bignum-ref result j)
1605                              (%ashl (%make-ones (- end-bits pos-bits))
1606                                     pos-bits))))))
1607           (bignum-plusp)
1608           (t
1609            (setf (%bignum-ref result j)
1610                  (%ashl (%make-ones bits-past-pos-bits) pos-bits))
1611            (do ((j j+1 (1+ j)))
1612                ((= j end-digit)
1613                 (setf (%bignum-ref result j) (%make-ones end-bits)))
1614              (declare (type bignum-index j))
1615              (setf (%bignum-ref result j) all-ones-digit)))))
1616       (declare (type bignum-index i j j+1))
1617       (let ((digit (%bignum-ref bignum-byte i)))
1618         (declare (type bignum-element-type digit))
1619         (setf (%bignum-ref result j)
1620               (logior (%bignum-ref result j)
1621                       (%ashl (logand digit high-mask) pos-bits)))
1622         (setf (%bignum-ref result j+1)
1623               (logand (ash digit minus-high-bits)
1624                       ;; LOGAND should be unnecessary here with a logical right
1625                       ;; shift or a correct unsigned-byte-32 one.
1626                       low-mask))))))
1627 |#
1628 \f
1629 ;;;; TRUNCATE
1630
1631 ;;; This is the original sketch of the algorithm from which I implemented this
1632 ;;; TRUNCATE, assuming both operands are bignums. I should modify this to work
1633 ;;; with the documentation on my functions, as a general introduction. I've
1634 ;;; left this here just in case someone needs it in the future. Don't look at
1635 ;;; this unless reading the functions' comments leaves you at a loss. Remember
1636 ;;; this comes from Knuth, so the book might give you the right general
1637 ;;; overview.
1638 ;;;
1639 ;;; (truncate x y):
1640 ;;;
1641 ;;; If X's magnitude is less than Y's, then result is 0 with remainder X.
1642 ;;;
1643 ;;; Make x and y positive, copying x if it is already positive.
1644 ;;;
1645 ;;; Shift y left until there's a 1 in the 30'th bit (most significant, non-sign
1646 ;;;       digit)
1647 ;;;    Just do most sig digit to determine how much to shift whole number.
1648 ;;; Shift x this much too.
1649 ;;; Remember this initial shift count.
1650 ;;;
1651 ;;; Allocate q to be len-x minus len-y quantity plus 1.
1652 ;;;
1653 ;;; i = last digit of x.
1654 ;;; k = last digit of q.
1655 ;;;
1656 ;;; LOOP
1657 ;;;
1658 ;;; j = last digit of y.
1659 ;;;
1660 ;;; compute guess.
1661 ;;; if x[i] = y[j] then g = #xFFFFFFFF
1662 ;;; else g = x[i]x[i-1]/y[j].
1663 ;;;
1664 ;;; check guess.
1665 ;;; %UNSIGNED-MULTIPLY returns b and c defined below.
1666 ;;;    a = x[i-1] - (logand (* g y[j]) #xFFFFFFFF).
1667 ;;;       Use %UNSIGNED-MULTIPLY taking low-order result.
1668 ;;;    b = (logand (ash (* g y[j-1]) -32) #xFFFFFFFF).
1669 ;;;    c = (logand (* g y[j-1]) #xFFFFFFFF).
1670 ;;; if a < b, okay.
1671 ;;; if a > b, guess is too high
1672 ;;;    g = g - 1; go back to "check guess".
1673 ;;; if a = b and c > x[i-2], guess is too high
1674 ;;;    g = g - 1; go back to "check guess".
1675 ;;; GUESS IS 32-BIT NUMBER, SO USE THING TO KEEP IN SPECIAL REGISTER
1676 ;;; SAME FOR A, B, AND C.
1677 ;;;
1678 ;;; Subtract g * y from x[i - len-y+1]..x[i]. See paper for doing this in step.
1679 ;;; If x[i] < 0, guess is screwed up.
1680 ;;;    negative g, then add 1
1681 ;;;    zero or positive g, then subtract 1
1682 ;;; AND add y back into x[len-y+1..i].
1683 ;;;
1684 ;;; q[k] = g.
1685 ;;; i = i - 1.
1686 ;;; k = k - 1.
1687 ;;;
1688 ;;; If k>=0, goto LOOP.
1689 ;;;
1690 ;;; Now quotient is good, but remainder is not.
1691 ;;; Shift x right by saved initial left shifting count.
1692 ;;;
1693 ;;; Check quotient and remainder signs.
1694 ;;; x pos y pos --> q pos r pos
1695 ;;; x pos y neg --> q neg r pos
1696 ;;; x neg y pos --> q neg r neg
1697 ;;; x neg y neg --> q pos r neg
1698 ;;;
1699 ;;; Normalize quotient and remainder. Cons result if necessary.
1700
1701 ;;; These are used by BIGNUM-TRUNCATE and friends in the general case.
1702 (defvar *truncate-x*)
1703 (defvar *truncate-y*)
1704
1705 ;;; This divides x by y returning the quotient and remainder. In the general
1706 ;;; case, we shift y to setup for the algorithm, and we use two buffers to save
1707 ;;; consing intermediate values. X gets destructively modified to become the
1708 ;;; remainder, and we have to shift it to account for the initial Y shift.
1709 ;;; After we multiple bind q and r, we first fix up the signs and then return
1710 ;;; the normalized results.
1711 (defun bignum-truncate (x y)
1712   (declare (type bignum-type x y))
1713   (let* ((x-plusp (%bignum-0-or-plusp x (%bignum-length x)))
1714          (y-plusp (%bignum-0-or-plusp y (%bignum-length y)))
1715          (x (if x-plusp x (negate-bignum x nil)))
1716          (y (if y-plusp y (negate-bignum y nil)))
1717          (len-x (%bignum-length x))
1718          (len-y (%bignum-length y)))
1719     (multiple-value-bind (q r)
1720         (cond ((< len-y 2)
1721                (bignum-truncate-single-digit x len-x y))
1722               ((plusp (bignum-compare y x))
1723                (let ((res (%allocate-bignum len-x)))
1724                  (dotimes (i len-x)
1725                    (setf (%bignum-ref res i) (%bignum-ref x i)))
1726                  (values 0 res)))
1727               (t
1728                (let ((len-x+1 (1+ len-x)))
1729                  (with-bignum-buffers ((*truncate-x* len-x+1)
1730                                        (*truncate-y* (1+ len-y)))
1731                    (let ((y-shift (shift-y-for-truncate y)))
1732                      (shift-and-store-truncate-buffers x len-x y len-y y-shift)
1733                      (values (do-truncate len-x+1 len-y)
1734                              ;; DO-TRUNCATE must execute first.
1735                              (cond
1736                               ((zerop y-shift)
1737                                (let ((res (%allocate-bignum len-y)))
1738                                  (declare (type bignum-type res))
1739                                  (bignum-replace res *truncate-x* :end2 len-y)
1740                                  (%normalize-bignum res len-y)))
1741                               (t
1742                                (shift-right-unaligned
1743                                 *truncate-x* 0 y-shift len-y
1744                                 ((= j res-len-1)
1745                                  (setf (%bignum-ref res j)
1746                                        (%ashr (%bignum-ref *truncate-x* i)
1747                                               y-shift))
1748                                  (%normalize-bignum res res-len))
1749                                 res)))))))))
1750       (let ((quotient (cond ((eq x-plusp y-plusp) q)
1751                             ((typep q 'fixnum) (the fixnum (- q)))
1752                             (t (negate-bignum-in-place q))))
1753             (rem (cond (x-plusp r)
1754                        ((typep r 'fixnum) (the fixnum (- r)))
1755                        (t (negate-bignum-in-place r)))))
1756         (values (if (typep quotient 'fixnum)
1757                     quotient
1758                     (%normalize-bignum quotient (%bignum-length quotient)))
1759                 (if (typep rem 'fixnum)
1760                     rem
1761                     (%normalize-bignum rem (%bignum-length rem))))))))
1762
1763 ;;; This divides x by y when y is a single bignum digit. BIGNUM-TRUNCATE fixes
1764 ;;; up the quotient and remainder with respect to sign and normalization.
1765 ;;;
1766 ;;; We don't have to worry about shifting y to make its most significant digit
1767 ;;; sufficiently large for %FLOOR to return 32-bit quantities for the q-digit
1768 ;;; and r-digit. If y is a single digit bignum, it is already large enough
1769 ;;; for %FLOOR. That is, it has some bits on pretty high in the digit.
1770 (defun bignum-truncate-single-digit (x len-x y)
1771   (declare (type bignum-index len-x))
1772   (let ((q (%allocate-bignum len-x))
1773         (r 0)
1774         (y (%bignum-ref y 0)))
1775     (declare (type bignum-element-type r y))
1776     (do ((i (1- len-x) (1- i)))
1777         ((minusp i))
1778       (multiple-value-bind (q-digit r-digit) (%floor r (%bignum-ref x i) y)
1779         (declare (type bignum-element-type q-digit r-digit))
1780         (setf (%bignum-ref q i) q-digit)
1781         (setf r r-digit)))
1782     (let ((rem (%allocate-bignum 1)))
1783       (setf (%bignum-ref rem 0) r)
1784       (values q rem))))
1785
1786 ;;; This divides *truncate-x* by *truncate-y*, and len-x and len-y tell us how
1787 ;;; much of the buffers we care about. TRY-BIGNUM-TRUNCATE-GUESS modifies
1788 ;;; *truncate-x* on each interation, and this buffer becomes our remainder.
1789 ;;;
1790 ;;; *truncate-x* definitely has at least three digits, and it has one more than
1791 ;;; *truncate-y*. This keeps i, i-1, i-2, and low-x-digit happy. Thanks to
1792 ;;; SHIFT-AND-STORE-TRUNCATE-BUFFERS.
1793 (defun do-truncate (len-x len-y)
1794   (declare (type bignum-index len-x len-y))
1795   (let* ((len-q (- len-x len-y))
1796          ;; Add one for extra sign digit in case high bit is on.
1797          (q (%allocate-bignum (1+ len-q)))
1798          (k (1- len-q))
1799          (y1 (%bignum-ref *truncate-y* (1- len-y)))
1800          (y2 (%bignum-ref *truncate-y* (- len-y 2)))
1801          (i (1- len-x))
1802          (i-1 (1- i))
1803          (i-2 (1- i-1))
1804          (low-x-digit (- i len-y)))
1805     (declare (type bignum-index len-q k i i-1 i-2 low-x-digit)
1806              (type bignum-element-type y1 y2))
1807     (loop
1808       (setf (%bignum-ref q k)
1809             (try-bignum-truncate-guess
1810              ;; This modifies *truncate-x*. Must access elements each pass.
1811              (bignum-truncate-guess y1 y2
1812                                     (%bignum-ref *truncate-x* i)
1813                                     (%bignum-ref *truncate-x* i-1)
1814                                     (%bignum-ref *truncate-x* i-2))
1815              len-y low-x-digit))
1816       (cond ((zerop k) (return))
1817             (t (decf k)
1818                (decf low-x-digit)
1819                (shiftf i i-1 i-2 (1- i-2)))))
1820     q))
1821
1822 ;;; This takes a digit guess, multiplies it by *truncate-y* for a result one
1823 ;;; greater in length than len-y, and subtracts this result from *truncate-x*.
1824 ;;; Low-x-digit is the first digit of x to start the subtraction, and we know x
1825 ;;; is long enough to subtract a len-y plus one length bignum from it. Next we
1826 ;;; check the result of the subtraction, and if the high digit in x became
1827 ;;; negative, then our guess was one too big. In this case, return one less
1828 ;;; than guess passed in, and add one value of y back into x to account for
1829 ;;; subtracting one too many. Knuth shows that the guess is wrong on the order
1830 ;;; of 3/b, where b is the base (2 to the digit-size power) -- pretty rarely.
1831 (defun try-bignum-truncate-guess (guess len-y low-x-digit)
1832   (declare (type bignum-index low-x-digit len-y)
1833            (type bignum-element-type guess))
1834   (let ((carry-digit 0)
1835         (borrow 1)
1836         (i low-x-digit))
1837     (declare (type bignum-element-type carry-digit)
1838              (type bignum-index i)
1839              (fixnum borrow))
1840     ;; Multiply guess and divisor, subtracting from dividend simultaneously.
1841     (dotimes (j len-y)
1842       (multiple-value-bind (high-digit low-digit)
1843           (%multiply-and-add guess
1844                              (%bignum-ref *truncate-y* j)
1845                              carry-digit)
1846         (declare (type bignum-element-type high-digit low-digit))
1847         (setf carry-digit high-digit)
1848         (multiple-value-bind (x temp-borrow)
1849             (%subtract-with-borrow (%bignum-ref *truncate-x* i)
1850                                    low-digit
1851                                    borrow)
1852           (declare (type bignum-element-type x)
1853                    (fixnum temp-borrow))
1854           (setf (%bignum-ref *truncate-x* i) x)
1855           (setf borrow temp-borrow)))
1856       (incf i))
1857     (setf (%bignum-ref *truncate-x* i)
1858           (%subtract-with-borrow (%bignum-ref *truncate-x* i)
1859                                  carry-digit borrow))
1860     ;; See whether guess is off by one, adding one Y back in if necessary.
1861     (cond ((%digit-0-or-plusp (%bignum-ref *truncate-x* i))
1862            guess)
1863           (t
1864            ;; If subtraction has negative result, add one divisor value back
1865            ;; in. The guess was one too large in magnitude.
1866            (let ((i low-x-digit)
1867                  (carry 0))
1868              (dotimes (j len-y)
1869                (multiple-value-bind (v k)
1870                    (%add-with-carry (%bignum-ref *truncate-y* j)
1871                                     (%bignum-ref *truncate-x* i)
1872                                     carry)
1873                  (declare (type bignum-element-type v))
1874                  (setf (%bignum-ref *truncate-x* i) v)
1875                  (setf carry k))
1876                (incf i))
1877              (setf (%bignum-ref *truncate-x* i)
1878                    (%add-with-carry (%bignum-ref *truncate-x* i) 0 carry)))
1879            (%subtract-with-borrow guess 1 1)))))
1880
1881 ;;; This returns a guess for the next division step. Y1 is the highest y
1882 ;;; digit, and y2 is the second to highest y digit. The x... variables are
1883 ;;; the three highest x digits for the next division step.
1884 ;;;
1885 ;;; From Knuth, our guess is either all ones or x-i and x-i-1 divided by y1,
1886 ;;; depending on whether x-i and y1 are the same. We test this guess by
1887 ;;; determining whether guess*y2 is greater than the three high digits of x
1888 ;;; minus guess*y1 shifted left one digit:
1889 ;;;    ------------------------------
1890 ;;;   |    x-i    |   x-i-1  | x-i-2 |
1891 ;;;    ------------------------------
1892 ;;;    ------------------------------
1893 ;;; - | g*y1 high | g*y1 low |   0   |
1894 ;;;    ------------------------------
1895 ;;;             ...               <   guess*y2     ???
1896 ;;; If guess*y2 is greater, then we decrement our guess by one and try again.
1897 ;;; This returns a guess that is either correct or one too large.
1898 (defun bignum-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
1899   (declare (type bignum-element-type y1 y2 x-i x-i-1 x-i-2))
1900   (let ((guess (if (%digit-compare x-i y1)
1901                    all-ones-digit
1902                    (%floor x-i x-i-1 y1))))
1903     (declare (type bignum-element-type guess))
1904     (loop
1905       (multiple-value-bind (high-guess*y1 low-guess*y1) (%multiply guess y1)
1906         (declare (type bignum-element-type low-guess*y1 high-guess*y1))
1907         (multiple-value-bind (high-guess*y2 low-guess*y2)
1908             (%multiply guess y2)
1909           (declare (type bignum-element-type high-guess*y2 low-guess*y2))
1910           (multiple-value-bind (middle-digit borrow)
1911               (%subtract-with-borrow x-i-1 low-guess*y1 1)
1912             (declare (type bignum-element-type middle-digit)
1913                      (fixnum borrow))
1914             ;; Supplying borrow of 1 means there was no borrow, and we know
1915             ;; x-i-2 minus 0 requires no borrow.
1916             (let ((high-digit (%subtract-with-borrow x-i high-guess*y1 borrow)))
1917               (declare (type bignum-element-type high-digit))
1918               (if (and (%digit-compare high-digit 0)
1919                        (or (%digit-greater high-guess*y2 middle-digit)
1920                            (and (%digit-compare middle-digit high-guess*y2)
1921                                 (%digit-greater low-guess*y2 x-i-2))))
1922                   (setf guess (%subtract-with-borrow guess 1 1))
1923                   (return guess)))))))))
1924
1925 ;;; This returns the amount to shift y to place a one in the second highest
1926 ;;; bit. Y must be positive. If the last digit of y is zero, then y has a
1927 ;;; one in the previous digit's sign bit, so we know it will take one less
1928 ;;; than digit-size to get a one where we want. Otherwise, we count how many
1929 ;;; right shifts it takes to get zero; subtracting this value from digit-size
1930 ;;; tells us how many high zeros there are which is one more than the shift
1931 ;;; amount sought.
1932 ;;;
1933 ;;; Note: This is exactly the same as one less than the integer-length of the
1934 ;;; last digit subtracted from the digit-size.
1935 ;;;
1936 ;;; We shift y to make it sufficiently large that doing the 64-bit by 32-bit
1937 ;;; %FLOOR calls ensures the quotient and remainder fit in 32-bits.
1938 (defun shift-y-for-truncate (y)
1939   (let* ((len (%bignum-length y))
1940          (last (%bignum-ref y (1- len))))
1941     (declare (type bignum-index len)
1942              (type bignum-element-type last))
1943     (- digit-size (integer-length last) 1)))
1944
1945 ;;; Stores two bignums into the truncation bignum buffers, shifting them on the
1946 ;;; way in. This assumes x and y are positive and at least two in length, and
1947 ;;; it assumes *truncate-x* and *truncate-y* are one digit longer than x and y.
1948 (defun shift-and-store-truncate-buffers (x len-x y len-y shift)
1949   (declare (type bignum-index len-x len-y)
1950            (type (integer 0 (#.digit-size)) shift))
1951   (cond ((zerop shift)
1952          (bignum-replace *truncate-x* x :end1 len-x)
1953          (bignum-replace *truncate-y* y :end1 len-y))
1954         (t
1955          (bignum-ashift-left-unaligned x 0 shift (1+ len-x) *truncate-x*)
1956          (bignum-ashift-left-unaligned y 0 shift (1+ len-y) *truncate-y*))))
1957 \f
1958 ;;;; %FLOOR primitive for BIGNUM-TRUNCATE
1959
1960 ;;; When a machine leaves out a 64-bit by 32-bit divide instruction (that is,
1961 ;;; two bignum-digits divided by one), we have to roll our own (the hard way).
1962 ;;; Basically, we treat the operation as four 16-bit digits divided by two
1963 ;;; 16-bit digits. This means we have duplicated most of the code above to do
1964 ;;; this nearly general 16-bit digit bignum divide, but we've unrolled loops
1965 ;;; and made use of other properties of this specific divide situation.
1966
1967 ;;;; %FLOOR for machines with a 32x32 divider.
1968
1969 #!-sb-fluid
1970 (declaim (inline 32x16-subtract-with-borrow 32x16-add-with-carry
1971                  32x16-divide 32x16-multiply 32x16-multiply-split))
1972
1973 #!+32x16-divide
1974 (defconstant 32x16-base-1 #xFFFF)
1975
1976 ;;; This is similar to %SUBTRACT-WITH-BORROW. It returns a 16-bit difference
1977 ;;; and a borrow. Returning a 1 for the borrow means there was no borrow, and
1978 ;;; 0 means there was one.
1979 #!+32x16-divide
1980 (defun 32x16-subtract-with-borrow (a b borrow)
1981   (declare (type (unsigned-byte 16) a b)
1982            (type (integer 0 1) borrow))
1983   (let ((diff (+ (- a b) borrow 32x16-base-1)))
1984     (declare (type (unsigned-byte 17) diff))
1985     (values (logand diff #xFFFF)
1986             (ash diff -16))))
1987
1988 ;;; This adds a and b, 16-bit quantities, with the carry k. It returns a
1989 ;;; 16-bit sum and a second value, 0 or 1, indicating whether there was a
1990 ;;; carry.
1991 #!+32x16-divide
1992 (defun 32x16-add-with-carry (a b k)
1993   (declare (type (unsigned-byte 16) a b)
1994            (type (integer 0 1) k))
1995   (let ((res (the fixnum (+ a b k))))
1996     (declare (type (unsigned-byte 17) res))
1997     (if (zerop (the fixnum (logand #x10000 res)))
1998         (values res 0)
1999         (values (the (unsigned-byte 16) (logand #xFFFF res))
2000                 1))))
2001
2002 ;;; This is probably a 32-bit by 32-bit divide instruction.
2003 #!+32x16-divide
2004 (defun 32x16-divide (a b c)
2005   (declare (type (unsigned-byte 16) a b c))
2006   (floor (the bignum-element-type
2007               (logior (the bignum-element-type (ash a 16))
2008                       b))
2009          c))
2010
2011 ;;; This basically exists since we know the answer won't overflow
2012 ;;; bignum-element-type. It's probably just a basic multiply instruction, but
2013 ;;; it can't cons an intermediate bignum. The result goes in a non-descriptor
2014 ;;; register.
2015 #!+32x16-divide
2016 (defun 32x16-multiply (a b)
2017   (declare (type (unsigned-byte 16) a b))
2018   (the bignum-element-type (* a b)))
2019
2020 ;;; This multiplies a and b, 16-bit quantities, and returns the result as two
2021 ;;; 16-bit quantities, high and low.
2022 #!+32x16-divide
2023 (defun 32x16-multiply-split (a b)
2024   (let ((res (32x16-multiply a b)))
2025     (declare (the bignum-element-type res))
2026     (values (the (unsigned-byte 16) (logand #xFFFF (ash res -16)))
2027             (the (unsigned-byte 16) (logand #xFFFF res)))))
2028
2029 ;;; The %FLOOR below uses this buffer the same way BIGNUM-TRUNCATE uses
2030 ;;; *truncate-x*. There's no y buffer since we pass around the two 16-bit
2031 ;;; digits and use them slightly differently than the general truncation
2032 ;;; algorithm above.
2033 #!+32x16-divide
2034 (defvar *32x16-truncate-x* (make-array 4 :element-type '(unsigned-byte 16)
2035                                        :initial-element 0))
2036
2037 ;;; This does the same thing as the %FLOOR above, but it does it at Lisp level
2038 ;;; when there is no 64x32-bit divide instruction on the machine.
2039 ;;;
2040 ;;; It implements the higher level tactics of BIGNUM-TRUNCATE, but it makes use
2041 ;;; of special situation provided, four 16-bit digits divided by two 16-bit
2042 ;;; digits.
2043 #!+32x16-divide
2044 (defun %floor (a b c)
2045   (declare (type bignum-element-type a b c))
2046   ;; Setup *32x16-truncate-x* buffer from a and b.
2047   (setf (aref *32x16-truncate-x* 0)
2048         (the (unsigned-byte 16) (logand #xFFFF b)))
2049   (setf (aref *32x16-truncate-x* 1)
2050         (the (unsigned-byte 16)
2051              (logand #xFFFF
2052                      (the (unsigned-byte 16) (ash b -16)))))
2053   (setf (aref *32x16-truncate-x* 2)
2054         (the (unsigned-byte 16) (logand #xFFFF a)))
2055   (setf (aref *32x16-truncate-x* 3)
2056         (the (unsigned-byte 16)
2057              (logand #xFFFF
2058                      (the (unsigned-byte 16) (ash a -16)))))
2059   ;; From DO-TRUNCATE, but unroll the loop.
2060   (let* ((y1 (logand #xFFFF (ash c -16)))
2061          (y2 (logand #xFFFF c))
2062          (q (the bignum-element-type
2063                  (ash (32x16-try-bignum-truncate-guess
2064                        (32x16-truncate-guess y1 y2
2065                                              (aref *32x16-truncate-x* 3)
2066                                              (aref *32x16-truncate-x* 2)
2067                                              (aref *32x16-truncate-x* 1))
2068                        y1 y2 1)
2069                       16))))
2070     (declare (type bignum-element-type q)
2071              (type (unsigned-byte 16) y1 y2))
2072     (values (the bignum-element-type
2073                  (logior q
2074                          (the (unsigned-byte 16)
2075                               (32x16-try-bignum-truncate-guess
2076                                (32x16-truncate-guess
2077                                 y1 y2
2078                                 (aref *32x16-truncate-x* 2)
2079                                 (aref *32x16-truncate-x* 1)
2080                                 (aref *32x16-truncate-x* 0))
2081                                y1 y2 0))))
2082             (the bignum-element-type
2083                  (logior (the bignum-element-type
2084                               (ash (aref *32x16-truncate-x* 1) 16))
2085                          (the (unsigned-byte 16)
2086                               (aref *32x16-truncate-x* 0)))))))
2087
2088 ;;; This is similar to TRY-BIGNUM-TRUNCATE-GUESS, but this unrolls the two
2089 ;;; loops. This also substitutes for %DIGIT-0-OR-PLUSP the equivalent
2090 ;;; expression without any embellishment or pretense of abstraction. The first
2091 ;;; loop is unrolled, but we've put the body of the loop into the function
2092 ;;; 32X16-TRY-GUESS-ONE-RESULT-DIGIT.
2093 #!+32x16-divide
2094 (defun 32x16-try-bignum-truncate-guess (guess y-high y-low low-x-digit)
2095   (declare (type bignum-index low-x-digit)
2096            (type (unsigned-byte 16) guess y-high y-low))
2097   (let ((high-x-digit (+ 2 low-x-digit)))
2098     ;; Multiply guess and divisor, subtracting from dividend simultaneously.
2099     (multiple-value-bind (guess*y-hold carry borrow)
2100         (32x16-try-guess-one-result-digit guess y-low 0 0 1 low-x-digit)
2101       (declare (type (unsigned-byte 16) guess*y-hold)
2102                (fixnum carry borrow))
2103       (multiple-value-bind (guess*y-hold carry borrow)
2104           (32x16-try-guess-one-result-digit guess y-high guess*y-hold
2105                                             carry borrow (1+ low-x-digit))
2106         (declare (type (unsigned-byte 16) guess*y-hold)
2107                  (fixnum borrow)
2108                  (ignore carry))
2109         (setf (aref *32x16-truncate-x* high-x-digit)
2110               (32x16-subtract-with-borrow (aref *32x16-truncate-x* high-x-digit)
2111                                           guess*y-hold borrow))))
2112     ;; See whether guess is off by one, adding one Y back in if necessary.
2113     (cond ((zerop (logand #x8000 (aref *32x16-truncate-x* high-x-digit)))
2114            ;; The subtraction result is zero or positive.
2115            guess)
2116           (t
2117            ;; If subtraction has negative result, add one divisor value back
2118            ;; in. The guess was one too large in magnitude.
2119            (multiple-value-bind (v carry)
2120                (32x16-add-with-carry y-low
2121                                      (aref *32x16-truncate-x* low-x-digit)
2122                                      0)
2123              (declare (type (unsigned-byte 16) v))
2124              (setf (aref *32x16-truncate-x* low-x-digit) v)
2125              (multiple-value-bind (v carry)
2126                  (32x16-add-with-carry y-high
2127                                        (aref *32x16-truncate-x*
2128                                              (1+ low-x-digit))
2129                                        carry)
2130                (setf (aref *32x16-truncate-x* (1+ low-x-digit)) v)
2131                (setf (aref *32x16-truncate-x* high-x-digit)
2132                      (32x16-add-with-carry (aref *32x16-truncate-x* high-x-digit)
2133                                            carry 0))))
2134            (if (zerop (logand #x8000 guess))
2135                (1- guess)
2136                (1+ guess))))))
2137
2138 ;;; This is similar to the body of the loop in TRY-BIGNUM-TRUNCATE-GUESS that
2139 ;;; multiplies the guess by y and subtracts the result from x simultaneously.
2140 ;;; This returns the digit remembered as part of the multiplication, the carry
2141 ;;; from additions done on behalf of the multiplication, and the borrow from
2142 ;;; doing the subtraction.
2143 #!+32x16-divide
2144 (defun 32x16-try-guess-one-result-digit (guess y-digit guess*y-hold
2145                                          carry borrow x-index)
2146   (multiple-value-bind (high-digit low-digit)
2147       (32x16-multiply-split guess y-digit)
2148     (declare (type (unsigned-byte 16) high-digit low-digit))
2149     (multiple-value-bind (low-digit temp-carry)
2150         (32x16-add-with-carry low-digit guess*y-hold carry)
2151       (declare (type (unsigned-byte 16) low-digit))
2152       (multiple-value-bind (high-digit temp-carry)
2153           (32x16-add-with-carry high-digit temp-carry 0)
2154         (declare (type (unsigned-byte 16) high-digit))
2155         (multiple-value-bind (x temp-borrow)
2156             (32x16-subtract-with-borrow (aref *32x16-truncate-x* x-index)
2157                                         low-digit borrow)
2158           (declare (type (unsigned-byte 16) x))
2159           (setf (aref *32x16-truncate-x* x-index) x)
2160           (values high-digit temp-carry temp-borrow))))))
2161
2162 ;;; This is similar to BIGNUM-TRUNCATE-GUESS, but instead of computing the
2163 ;;; guess exactly as described in the its comments (digit by digit), this
2164 ;;; massages the 16-bit quantities into 32-bit quantities and performs the
2165 #!+32x16-divide
2166 (defun 32x16-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
2167   (declare (type (unsigned-byte 16) y1 y2 x-i x-i-1 x-i-2))
2168   (let ((guess (if (= x-i y1)
2169                    #xFFFF
2170                    (32x16-divide x-i x-i-1 y1))))
2171     (declare (type (unsigned-byte 16) guess))
2172     (loop
2173       (let* ((guess*y1 (the bignum-element-type
2174                             (ash (logand #xFFFF
2175                                          (the bignum-element-type
2176                                               (32x16-multiply guess y1)))
2177                                  16)))
2178              (x-y (%subtract-with-borrow
2179                    (the bignum-element-type
2180                         (logior (the bignum-element-type
2181                                      (ash x-i-1 16))
2182                                 x-i-2))
2183                    guess*y1
2184                    1))
2185              (guess*y2 (the bignum-element-type (%multiply guess y2))))
2186         (declare (type bignum-element-type guess*y1 x-y guess*y2))
2187         (if (%digit-greater guess*y2 x-y)
2188             (decf guess)
2189             (return guess))))))
2190 \f
2191 ;;;; general utilities
2192
2193 ;;; Allocate a single word bignum that holds fixnum. This is useful when
2194 ;;; we are trying to mix fixnum and bignum operands.
2195 #!-sb-fluid (declaim (inline make-small-bignum))
2196 (defun make-small-bignum (fixnum)
2197   (let ((res (%allocate-bignum 1)))
2198     (setf (%bignum-ref res 0) (%fixnum-to-digit fixnum))
2199     res))
2200
2201 ;;; Internal in-place operations use this to fixup remaining digits in the
2202 ;;; incoming data, such as in-place shifting. This is basically the same as
2203 ;;; the first form in %NORMALIZE-BIGNUM, but we return the length of the buffer
2204 ;;; instead of shrinking the bignum.
2205 #!-sb-fluid (declaim (sb!ext:maybe-inline %normalize-bignum-buffer))
2206 (defun %normalize-bignum-buffer (result len)
2207   (declare (type bignum-type result)
2208            (type bignum-index len))
2209   (unless (= len 1)
2210     (do ((next-digit (%bignum-ref result (- len 2))
2211                      (%bignum-ref result (- len 2)))
2212          (sign-digit (%bignum-ref result (1- len)) next-digit))
2213         ((not (zerop (logxor sign-digit (%ashr next-digit (1- digit-size))))))
2214         (decf len)
2215         (setf (%bignum-ref result len) 0)       
2216         (when (= len 1)
2217               (return))))
2218   len)
2219
2220 ;;; This drops the last digit if it is unnecessary sign information. It repeats
2221 ;;; this as needed, possibly ending with a fixnum. If the resulting length from
2222 ;;; shrinking is one, see whether our one word is a fixnum. Shift the possible
2223 ;;; fixnum bits completely out of the word, and compare this with shifting the
2224 ;;; sign bit all the way through. If the bits are all 1's or 0's in both words,
2225 ;;; then there are just sign bits between the fixnum bits and the sign bit. If
2226 ;;; we do have a fixnum, shift it over for the two low-tag bits.
2227 (defun %normalize-bignum (result len)
2228   (declare (type bignum-type result)
2229            (type bignum-index len)
2230            (inline %normalize-bignum-buffer))
2231   (let ((newlen (%normalize-bignum-buffer result len)))
2232     (declare (type bignum-index newlen))
2233     (unless (= newlen len)
2234       (%bignum-set-length result newlen))
2235     (if (= newlen 1)
2236         (let ((digit (%bignum-ref result 0)))
2237           (if (= (%ashr digit 29) (%ashr digit (1- digit-size)))
2238               (%fixnum-digit-with-correct-sign digit)
2239               result))
2240         result)))
2241
2242 ;;; This drops the last digit if it is unnecessary sign information. It
2243 ;;; repeats this as needed, possibly ending with a fixnum magnitude but never
2244 ;;; returning a fixnum.
2245 (defun %mostly-normalize-bignum (result len)
2246   (declare (type bignum-type result)
2247            (type bignum-index len)
2248            (inline %normalize-bignum-buffer))
2249   (let ((newlen (%normalize-bignum-buffer result len)))
2250     (declare (type bignum-index newlen))
2251     (unless (= newlen len)
2252       (%bignum-set-length result newlen))
2253     result))
2254 \f
2255 ;;;; hashing
2256
2257 ;;; the bignum case of the SXHASH function
2258 (defun sxhash-bignum (x)
2259   (let ((result 316495330))
2260     (declare (type fixnum result))
2261     (dotimes (i (%bignum-length x))
2262       (declare (type index i))
2263       (let ((xi (%bignum-ref x i)))
2264         (mixf result
2265               (logand most-positive-fixnum
2266                       xi
2267                       (ash xi -7)))))
2268     result))