1.0.11.34: better SUBSEQ on lists
[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-logbitp bignum-logcount
26 ;;;   These symbols define the interface to the compiler:
27 ;;;       bignum-type bignum-element-type bignum-index %allocate-bignum
28 ;;;       %bignum-length %bignum-set-length %bignum-ref %bignum-set
29 ;;;       %digit-0-or-plusp %add-with-carry %subtract-with-borrow
30 ;;;       %multiply-and-add %multiply %lognot %logand %logior %logxor
31 ;;;       %fixnum-to-digit %floor %fixnum-digit-with-correct-sign %ashl
32 ;;;       %ashr %digit-logical-shift-right))
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 machine
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 sb!vm:n-word-bits)
112
113 (defconstant maximum-bignum-length (1- (ash 1 (- sb!vm:n-word-bits
114                                                  sb!vm:n-widetag-bits))))
115
116 (defconstant all-ones-digit (1- (ash 1 sb!vm:n-word-bits)))
117 \f
118 ;;;; internal inline routines
119
120 ;;; %ALLOCATE-BIGNUM must zero all elements.
121 (defun %allocate-bignum (length)
122   (declare (type bignum-index length))
123   (%allocate-bignum length))
124
125 ;;; Extract the length of the bignum.
126 (defun %bignum-length (bignum)
127   (declare (type bignum-type bignum))
128   (%bignum-length bignum))
129
130 ;;; %BIGNUM-REF needs to access bignums as obviously as possible, and it needs
131 ;;; to be able to return the digit somewhere no one looks for real objects.
132 (defun %bignum-ref (bignum i)
133   (declare (type bignum-type bignum)
134            (type bignum-index i))
135   (%bignum-ref bignum i))
136 (defun %bignum-set (bignum i value)
137   (declare (type bignum-type bignum)
138            (type bignum-index i)
139            (type bignum-element-type value))
140   (%bignum-set bignum i value))
141
142 ;;; Return T if digit is positive, or NIL if negative.
143 (defun %digit-0-or-plusp (digit)
144   (declare (type bignum-element-type digit))
145   (not (logbitp (1- digit-size) digit)))
146
147 #!-sb-fluid (declaim (inline %bignum-0-or-plusp))
148 (defun %bignum-0-or-plusp (bignum len)
149   (declare (type bignum-type bignum)
150            (type bignum-index len))
151   (%digit-0-or-plusp (%bignum-ref bignum (1- len))))
152
153 ;;; This should be in assembler, and should not cons intermediate
154 ;;; results. It returns a bignum digit and a carry resulting from adding
155 ;;; together a, b, and an incoming carry.
156 (defun %add-with-carry (a b carry)
157   (declare (type bignum-element-type a b)
158            (type (mod 2) carry))
159   (%add-with-carry a b carry))
160
161 ;;; This should be in assembler, and should not cons intermediate
162 ;;; results. It returns a bignum digit and a borrow resulting from
163 ;;; subtracting b from a, and subtracting a possible incoming borrow.
164 ;;;
165 ;;; We really do:  a - b - 1 + borrow, where borrow is either 0 or 1.
166 (defun %subtract-with-borrow (a b borrow)
167   (declare (type bignum-element-type a b)
168            (type (mod 2) borrow))
169   (%subtract-with-borrow a b borrow))
170
171 ;;; Multiply two digit-size numbers, returning a 2*digit-size result
172 ;;; split into two digit-size quantities.
173 (defun %multiply (x y)
174   (declare (type bignum-element-type x y))
175   (%multiply x y))
176
177 ;;; This multiplies x-digit and y-digit, producing high and low digits
178 ;;; manifesting the result. Then it adds the low digit, res-digit, and
179 ;;; carry-in-digit. Any carries (note, you still have to add two digits
180 ;;; at a time possibly producing two carries) from adding these three
181 ;;; digits get added to the high digit from the multiply, producing the
182 ;;; next carry digit.  Res-digit is optional since two uses of this
183 ;;; primitive multiplies a single digit bignum by a multiple digit
184 ;;; bignum, and in this situation there is no need for a result buffer
185 ;;; accumulating partial results which is where the res-digit comes
186 ;;; from.
187 (defun %multiply-and-add (x-digit y-digit carry-in-digit
188                           &optional (res-digit 0))
189   (declare (type bignum-element-type x-digit y-digit res-digit carry-in-digit))
190   (%multiply-and-add x-digit y-digit carry-in-digit res-digit))
191
192 (defun %lognot (digit)
193   (declare (type bignum-element-type digit))
194   (%lognot digit))
195
196 ;;; Each of these does the digit-size unsigned op.
197 #!-sb-fluid (declaim (inline %logand %logior %logxor))
198 (defun %logand (a b)
199   (declare (type bignum-element-type a b))
200   (logand a b))
201 (defun %logior (a b)
202   (declare (type bignum-element-type a b))
203   (logior a b))
204 (defun %logxor (a b)
205   (declare (type bignum-element-type a b))
206   (logxor a b))
207
208 ;;; This takes a fixnum and sets it up as an unsigned digit-size
209 ;;; quantity.
210 (defun %fixnum-to-digit (x)
211   (declare (fixnum x))
212   (logand x (1- (ash 1 digit-size))))
213
214 #!-32x16-divide
215 ;;; This takes three digits and returns the FLOOR'ed result of
216 ;;; dividing the first two as a 2*digit-size integer by the third.
217 ;;;
218 ;;; Do weird LET and SETQ stuff to bamboozle the compiler into allowing
219 ;;; the %FLOOR transform to expand into pseudo-assembler for which the
220 ;;; compiler can later correctly allocate registers.
221 (defun %floor (a b c)
222   (let ((a a) (b b) (c c))
223     (declare (type bignum-element-type a b c))
224     (setq a a b b c c)
225     (%floor a b c)))
226
227 ;;; Convert the digit to a regular integer assuming that the digit is signed.
228 (defun %fixnum-digit-with-correct-sign (digit)
229   (declare (type bignum-element-type digit))
230   (if (logbitp (1- digit-size) digit)
231       (logior digit (ash -1 digit-size))
232       digit))
233
234 ;;; Do an arithmetic shift right of data even though bignum-element-type is
235 ;;; unsigned.
236 (defun %ashr (data count)
237   (declare (type bignum-element-type data)
238            (type (mod #.sb!vm:n-word-bits) count))
239   (%ashr data count))
240
241 ;;; This takes a digit-size quantity and shifts it to the left,
242 ;;; returning a digit-size quantity.
243 (defun %ashl (data count)
244   (declare (type bignum-element-type data)
245            (type (mod #.sb!vm:n-word-bits) count))
246   (%ashl data count))
247
248 ;;; Do an unsigned (logical) right shift of a digit by Count.
249 (defun %digit-logical-shift-right (data count)
250   (declare (type bignum-element-type data)
251            (type (mod #.sb!vm:n-word-bits) count))
252   (%digit-logical-shift-right data count))
253
254 ;;; Change the length of bignum to be newlen. Newlen must be the same or
255 ;;; smaller than the old length, and any elements beyond newlen must be zeroed.
256 (defun %bignum-set-length (bignum newlen)
257   (declare (type bignum-type bignum)
258            (type bignum-index newlen))
259   (%bignum-set-length bignum newlen))
260
261 ;;; This returns 0 or "-1" depending on whether the bignum is positive. This
262 ;;; is suitable for infinite sign extension to complete additions,
263 ;;; subtractions, negations, etc. This cannot return a -1 represented as
264 ;;; a negative fixnum since it would then have to low zeros.
265 #!-sb-fluid (declaim (inline %sign-digit))
266 (defun %sign-digit (bignum len)
267   (declare (type bignum-type bignum)
268            (type bignum-index len))
269   (%ashr (%bignum-ref bignum (1- len)) (1- digit-size)))
270
271 ;;; These take two digit-size quantities and compare or contrast them
272 ;;; without wasting time with incorrect type checking.
273 #!-sb-fluid (declaim (inline %digit-compare %digit-greater))
274 (defun %digit-compare (x y)
275   (= x y))
276 (defun %digit-greater (x y)
277   (> x y))
278 \f
279 (declaim (optimize (speed 3) (safety 0)))
280 \f
281 ;;;; addition
282
283 (defun add-bignums (a b)
284   (declare (type bignum-type a b))
285   (let ((len-a (%bignum-length a))
286         (len-b (%bignum-length b)))
287     (declare (type bignum-index len-a len-b))
288     (multiple-value-bind (a len-a b len-b)
289         (if (> len-a len-b)
290             (values a len-a b len-b)
291             (values b len-b a len-a))
292       (declare (type bignum-type a b)
293                (type bignum-index len-a len-b))
294       (let* ((len-res (1+ len-a))
295              (res (%allocate-bignum len-res))
296              (carry 0))
297         (declare (type bignum-index len-res)
298                  (type bignum-type res)
299                  (type (mod 2) carry))
300         (dotimes (i len-b)
301           (declare (type bignum-index i))
302           (multiple-value-bind (v k)
303               (%add-with-carry (%bignum-ref a i) (%bignum-ref b i) carry)
304             (declare (type bignum-element-type v)
305                      (type (mod 2) k))
306             (setf (%bignum-ref res i) v)
307             (setf carry k)))
308         (if (/= len-a len-b)
309             (finish-add a res carry (%sign-digit b len-b) len-b len-a)
310             (setf (%bignum-ref res len-a)
311                   (%add-with-carry (%sign-digit a len-a)
312                                    (%sign-digit b len-b)
313                                    carry)))
314         (%normalize-bignum res len-res)))))
315
316 ;;; This takes the longer of two bignums and propagates the carry through its
317 ;;; remaining high order digits.
318 (defun finish-add (a res carry sign-digit-b start end)
319   (declare (type bignum-type a res)
320            (type (mod 2) carry)
321            (type bignum-element-type sign-digit-b)
322            (type bignum-index start end))
323   (do ((i start (1+ i)))
324       ((= i end)
325        (setf (%bignum-ref res end)
326              (%add-with-carry (%sign-digit a end) sign-digit-b carry)))
327     (declare (type bignum-index i))
328     (multiple-value-bind (v k)
329         (%add-with-carry (%bignum-ref a i) sign-digit-b carry)
330       (setf (%bignum-ref res i) v)
331       (setf carry k)))
332   (values))
333 \f
334 ;;;; subtraction
335
336 (eval-when (:compile-toplevel :execute)
337
338 ;;; This subtracts b from a plugging result into res. Return-fun is the
339 ;;; function to call that fixes up the result returning any useful values, such
340 ;;; as the result. This macro may evaluate its arguments more than once.
341 (sb!xc:defmacro subtract-bignum-loop (a len-a b len-b res len-res return-fun)
342   (let ((borrow (gensym))
343         (a-digit (gensym))
344         (a-sign (gensym))
345         (b-digit (gensym))
346         (b-sign (gensym))
347         (i (gensym))
348         (v (gensym))
349         (k (gensym)))
350     `(let* ((,borrow 1)
351             (,a-sign (%sign-digit ,a ,len-a))
352             (,b-sign (%sign-digit ,b ,len-b)))
353        (declare (type bignum-element-type ,a-sign ,b-sign))
354        (dotimes (,i ,len-res)
355          (declare (type bignum-index ,i))
356          (let ((,a-digit (if (< ,i ,len-a) (%bignum-ref ,a ,i) ,a-sign))
357                (,b-digit (if (< ,i ,len-b) (%bignum-ref ,b ,i) ,b-sign)))
358            (declare (type bignum-element-type ,a-digit ,b-digit))
359            (multiple-value-bind (,v ,k)
360                (%subtract-with-borrow ,a-digit ,b-digit ,borrow)
361              (setf (%bignum-ref ,res ,i) ,v)
362              (setf ,borrow ,k))))
363        (,return-fun ,res ,len-res))))
364
365 ) ;EVAL-WHEN
366
367 (defun subtract-bignum (a b)
368   (declare (type bignum-type a b))
369   (let* ((len-a (%bignum-length a))
370          (len-b (%bignum-length b))
371          (len-res (1+ (max len-a len-b)))
372          (res (%allocate-bignum len-res)))
373     (declare (type bignum-index len-a len-b len-res)) ;Test len-res for bounds?
374     (subtract-bignum-loop a len-a b len-b res len-res %normalize-bignum)))
375
376 ;;; Operations requiring a subtraction without the overhead of intermediate
377 ;;; results, such as GCD, use this. It assumes Result is big enough for the
378 ;;; result.
379 (defun subtract-bignum-buffers-with-len (a len-a b len-b result len-res)
380   (declare (type bignum-type a b result)
381            (type bignum-index len-a len-b len-res))
382   (subtract-bignum-loop a len-a b len-b result len-res
383                         %normalize-bignum-buffer))
384
385 (defun subtract-bignum-buffers (a len-a b len-b result)
386   (declare (type bignum-type a b result)
387            (type bignum-index len-a len-b))
388   (subtract-bignum-loop a len-a b len-b result (max len-a len-b)
389                         %normalize-bignum-buffer))
390 \f
391 ;;;; multiplication
392
393 (defun multiply-bignums (a b)
394   (declare (type bignum-type a b))
395   (let* ((a-plusp (%bignum-0-or-plusp a (%bignum-length a)))
396          (b-plusp (%bignum-0-or-plusp b (%bignum-length b)))
397          (a (if a-plusp a (negate-bignum a)))
398          (b (if b-plusp b (negate-bignum b)))
399          (len-a (%bignum-length a))
400          (len-b (%bignum-length b))
401          (len-res (+ len-a len-b))
402          (res (%allocate-bignum len-res))
403          (negate-res (not (eq a-plusp b-plusp))))
404     (declare (type bignum-index len-a len-b len-res))
405     (dotimes (i len-a)
406       (declare (type bignum-index i))
407       (let ((carry-digit 0)
408             (x (%bignum-ref a i))
409             (k i))
410         (declare (type bignum-index k)
411                  (type bignum-element-type carry-digit x))
412         (dotimes (j len-b)
413           (multiple-value-bind (big-carry res-digit)
414               (%multiply-and-add x
415                                  (%bignum-ref b j)
416                                  (%bignum-ref res k)
417                                  carry-digit)
418             (declare (type bignum-element-type big-carry res-digit))
419             (setf (%bignum-ref res k) res-digit)
420             (setf carry-digit big-carry)
421             (incf k)))
422         (setf (%bignum-ref res k) carry-digit)))
423     (when negate-res (negate-bignum-in-place res))
424     (%normalize-bignum res len-res)))
425
426 (defun multiply-bignum-and-fixnum (bignum fixnum)
427   (declare (type bignum-type bignum) (type fixnum fixnum))
428   (let* ((bignum-plus-p (%bignum-0-or-plusp bignum (%bignum-length bignum)))
429          (fixnum-plus-p (not (minusp fixnum)))
430          (bignum (if bignum-plus-p bignum (negate-bignum bignum)))
431          (bignum-len (%bignum-length bignum))
432          (fixnum (if fixnum-plus-p fixnum (- fixnum)))
433          (result (%allocate-bignum (1+ bignum-len)))
434          (carry-digit 0))
435     (declare (type bignum-type bignum result)
436              (type bignum-index bignum-len)
437              (type bignum-element-type fixnum carry-digit))
438     (dotimes (index bignum-len)
439       (declare (type bignum-index index))
440       (multiple-value-bind (next-digit low)
441           (%multiply-and-add (%bignum-ref bignum index) fixnum carry-digit)
442         (declare (type bignum-element-type next-digit low))
443         (setf carry-digit next-digit)
444         (setf (%bignum-ref result index) low)))
445     (setf (%bignum-ref result bignum-len) carry-digit)
446     (unless (eq bignum-plus-p fixnum-plus-p)
447       (negate-bignum-in-place result))
448     (%normalize-bignum result (1+ bignum-len))))
449
450 (defun multiply-fixnums (a b)
451   (declare (fixnum a b))
452   (let* ((a-minusp (minusp a))
453          (b-minusp (minusp b)))
454     (multiple-value-bind (high low)
455         (%multiply (if a-minusp (- a) a)
456                    (if b-minusp (- b) b))
457       (declare (type bignum-element-type high low))
458       (if (and (zerop high)
459                (%digit-0-or-plusp low))
460           (let ((low (sb!ext:truly-the (unsigned-byte #.(1- sb!vm:n-word-bits))
461                                        (%fixnum-digit-with-correct-sign low))))
462             (if (eq a-minusp b-minusp)
463                 low
464                 (- low)))
465           (let ((res (%allocate-bignum 2)))
466             (%bignum-set res 0 low)
467             (%bignum-set res 1 high)
468             (unless (eq a-minusp b-minusp) (negate-bignum-in-place res))
469             (%normalize-bignum res 2))))))
470 \f
471 ;;;; BIGNUM-REPLACE and WITH-BIGNUM-BUFFERS
472
473 (eval-when (:compile-toplevel :execute)
474
475 (sb!xc:defmacro bignum-replace (dest
476                                 src
477                                 &key
478                                 (start1 '0)
479                                 end1
480                                 (start2 '0)
481                                 end2
482                                 from-end)
483   (sb!int:once-only ((n-dest dest)
484                      (n-src src))
485     (let ((n-start1 (gensym))
486           (n-end1 (gensym))
487           (n-start2 (gensym))
488           (n-end2 (gensym))
489           (i1 (gensym))
490           (i2 (gensym))
491           (end1 (or end1 `(%bignum-length ,n-dest)))
492           (end2 (or end2 `(%bignum-length ,n-src))))
493       (if from-end
494           `(let ((,n-start1 ,start1)
495                  (,n-start2 ,start2))
496              (do ((,i1 (1- ,end1) (1- ,i1))
497                   (,i2 (1- ,end2) (1- ,i2)))
498                  ((or (< ,i1 ,n-start1) (< ,i2 ,n-start2)))
499                (declare (fixnum ,i1 ,i2))
500                (%bignum-set ,n-dest ,i1
501                             (%bignum-ref ,n-src ,i2))))
502           (if (eql start1 start2)
503               `(let ((,n-end1 (min ,end1 ,end2)))
504                  (do ((,i1 ,start1 (1+ ,i1)))
505                      ((>= ,i1 ,n-end1))
506                    (declare (type bignum-index ,i1))
507                    (%bignum-set ,n-dest ,i1
508                                 (%bignum-ref ,n-src ,i1))))
509               `(let ((,n-end1 ,end1)
510                      (,n-end2 ,end2))
511                  (do ((,i1 ,start1 (1+ ,i1))
512                       (,i2 ,start2 (1+ ,i2)))
513                      ((or (>= ,i1 ,n-end1) (>= ,i2 ,n-end2)))
514                    (declare (type bignum-index ,i1 ,i2))
515                    (%bignum-set ,n-dest ,i1
516                                 (%bignum-ref ,n-src ,i2)))))))))
517
518 (sb!xc:defmacro with-bignum-buffers (specs &body body)
519   #!+sb-doc
520   "WITH-BIGNUM-BUFFERS ({(var size [init])}*) Form*"
521   (sb!int:collect ((binds)
522                    (inits))
523     (dolist (spec specs)
524       (let ((name (first spec))
525             (size (second spec)))
526         (binds `(,name (%allocate-bignum ,size)))
527         (let ((init (third spec)))
528           (when init
529             (inits `(bignum-replace ,name ,init))))))
530     `(let* ,(binds)
531        ,@(inits)
532        ,@body)))
533
534 ) ;EVAL-WHEN
535 \f
536 ;;;; GCD
537
538 (eval-when (:compile-toplevel :load-toplevel :execute)
539   ;; The asserts in the GCD implementation are way too expensive to
540   ;; check in normal use, and are disabled here.
541   (sb!xc:defmacro gcd-assert (&rest args)
542     (if nil
543         `(assert ,@args)))
544   ;; We'll be doing a lot of modular arithmetic.
545   (sb!xc:defmacro modularly (form)
546     `(logand all-ones-digit ,form)))
547
548 ;;; I'm not sure why I need this FTYPE declaration.  Compiled by the
549 ;;; target compiler, it can deduce the return type fine, but without
550 ;;; it, we pay a heavy price in BIGNUM-GCD when compiled by the
551 ;;; cross-compiler. -- CSR, 2004-07-19
552 (declaim (ftype (sfunction (bignum-type bignum-index bignum-type bignum-index)
553                            (and unsigned-byte fixnum))
554                 bignum-factors-of-two))
555 (defun bignum-factors-of-two (a len-a b len-b)
556   (declare (type bignum-index len-a len-b) (type bignum-type a b))
557   (do ((i 0 (1+ i))
558        (end (min len-a len-b)))
559       ((= i end) (error "Unexpected zero bignums?"))
560     (declare (type bignum-index i end))
561     (let ((or-digits (%logior (%bignum-ref a i) (%bignum-ref b i))))
562       (unless (zerop or-digits)
563         (return (do ((j 0 (1+ j))
564                      (or-digits or-digits (%ashr or-digits 1)))
565                     ((oddp or-digits) (+ (* i digit-size) j))
566                   (declare (type (mod #.sb!vm:n-word-bits) j))))))))
567
568 ;;; Multiply a bignum buffer with a fixnum or a digit, storing the
569 ;;; result in another bignum buffer, and without using any
570 ;;; temporaries. Inlined to avoid boxing smallnum if it's actually a
571 ;;; digit. Needed by GCD, should possibly OAOO with
572 ;;; MULTIPLY-BIGNUM-AND-FIXNUM.
573 (declaim (inline multiply-bignum-buffer-and-smallnum-to-buffer))
574 (defun multiply-bignum-buffer-and-smallnum-to-buffer (bignum bignum-len
575                                                              smallnum res)
576   (declare (type bignum-type bignum))
577   (let* ((bignum-plus-p (%bignum-0-or-plusp bignum bignum-len))
578          (smallnum-plus-p (not (minusp smallnum)))
579          (smallnum (if smallnum-plus-p smallnum (- smallnum)))
580          (carry-digit 0))
581     (declare (type bignum-type bignum res)
582              (type bignum-index bignum-len)
583              (type bignum-element-type smallnum carry-digit))
584     (unless bignum-plus-p
585       (negate-bignum-buffer-in-place bignum bignum-len))
586     (dotimes (index bignum-len)
587       (declare (type bignum-index index))
588       (multiple-value-bind (next-digit low)
589           (%multiply-and-add (%bignum-ref bignum index)
590                              smallnum
591                              carry-digit)
592         (declare (type bignum-element-type next-digit low))
593         (setf carry-digit next-digit)
594         (setf (%bignum-ref res index) low)))
595     (setf (%bignum-ref res bignum-len) carry-digit)
596     (unless bignum-plus-p
597       (negate-bignum-buffer-in-place bignum bignum-len))
598     (let ((res-len (%normalize-bignum-buffer res (1+ bignum-len))))
599       (unless (eq bignum-plus-p smallnum-plus-p)
600         (negate-bignum-buffer-in-place res res-len))
601       res-len)))
602
603 ;;; Given U and V, return U / V mod 2^32. Implements the algorithm in the
604 ;;; paper, but uses some clever bit-twiddling nicked from Nickle to do it.
605 (declaim (inline bmod))
606 (defun bmod (u v)
607   (let ((ud (%bignum-ref u 0))
608         (vd (%bignum-ref v 0))
609         (umask 0)
610         (imask 1)
611         (m 0))
612     (declare (type (unsigned-byte #.sb!vm:n-word-bits) ud vd umask imask m))
613     (dotimes (i digit-size)
614       (setf umask (logior umask imask))
615       (when (logtest ud umask)
616         (setf ud (modularly (- ud vd)))
617         (setf m (modularly (logior m imask))))
618       (setf imask (modularly (ash imask 1)))
619       (setf vd (modularly (ash vd 1))))
620     m))
621
622 (defun dmod (u u-len v v-len tmp1)
623   (loop while (> (bignum-buffer-integer-length u u-len)
624                  (+ (bignum-buffer-integer-length v v-len)
625                     digit-size))
626     do
627     (unless (zerop (%bignum-ref u 0))
628       (let* ((bmod (bmod u v))
629              (tmp1-len (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
630                                                                       bmod
631                                                                       tmp1)))
632         (setf u-len (subtract-bignum-buffers u u-len
633                                              tmp1 tmp1-len
634                                              u))
635         (bignum-abs-buffer u u-len)))
636     (gcd-assert (zerop (%bignum-ref u 0)))
637     (setf u-len (bignum-buffer-ashift-right u u-len digit-size)))
638   (let* ((d (+ 1 (- (bignum-buffer-integer-length u u-len)
639                     (bignum-buffer-integer-length v v-len))))
640          (n (1- (ash 1 d))))
641     (declare (type (unsigned-byte #.(integer-length #.sb!vm:n-word-bits)) d)
642              (type (unsigned-byte #.sb!vm:n-word-bits) n))
643     (gcd-assert (>= d 0))
644     (when (logtest (%bignum-ref u 0) n)
645       (let ((tmp1-len
646              (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
647                                                             (logand n (bmod u
648                                                                             v))
649                                                             tmp1)))
650         (setf u-len (subtract-bignum-buffers u u-len
651                                              tmp1 tmp1-len
652                                              u))
653         (bignum-abs-buffer u u-len)))
654     u-len))
655
656 (defconstant lower-ones-digit (1- (ash 1 (truncate sb!vm:n-word-bits 2))))
657
658 ;;; Find D and N such that (LOGAND ALL-ONES-DIGIT (- (* D X) (* N Y))) is 0,
659 ;;; (< 0 N LOWER-ONES-DIGIT) and (< 0 (ABS D) LOWER-ONES-DIGIT).
660 (defun reduced-ratio-mod (x y)
661   (let* ((c (bmod x y))
662          (n1 c)
663          (d1 1)
664          (n2 (modularly (1+ (modularly (lognot n1)))))
665          (d2 (modularly -1)))
666     (declare (type (unsigned-byte #.sb!vm:n-word-bits) n1 d1 n2 d2))
667     (loop while (> n2 (expt 2 (truncate digit-size 2))) do
668           (loop for i of-type (mod #.sb!vm:n-word-bits)
669                 downfrom (- (integer-length n1) (integer-length n2))
670                 while (>= n1 n2) do
671                 (when (>= n1 (modularly (ash n2 i)))
672                   (psetf n1 (modularly (- n1 (modularly (ash n2 i))))
673                          d1 (modularly (- d1 (modularly (ash d2 i)))))))
674           (psetf n1 n2
675                  d1 d2
676                  n2 n1
677                  d2 d1))
678     (values n2 (if (>= d2 (expt 2 (1- digit-size)))
679                    (lognot (logand most-positive-fixnum (lognot d2)))
680                    (logand lower-ones-digit d2)))))
681
682
683 (defun copy-bignum (a &optional (len (%bignum-length a)))
684   (let ((b (%allocate-bignum len)))
685     (bignum-replace b a)
686     (%bignum-set-length b len)
687     b))
688
689 ;;; Allocate a single word bignum that holds fixnum. This is useful when
690 ;;; we are trying to mix fixnum and bignum operands.
691 #!-sb-fluid (declaim (inline make-small-bignum))
692 (defun make-small-bignum (fixnum)
693   (let ((res (%allocate-bignum 1)))
694     (setf (%bignum-ref res 0) (%fixnum-to-digit fixnum))
695     res))
696
697 ;; When the larger number is less than this many bignum digits long, revert
698 ;; to old algorithm.
699 (defparameter *accelerated-gcd-cutoff* 3)
700
701 ;;; Alternate between k-ary reduction with the help of
702 ;;; REDUCED-RATIO-MOD and digit modulus reduction via DMOD. Once the
703 ;;; arguments get small enough, drop through to BIGNUM-MOD-GCD (since
704 ;;; k-ary reduction can introduce spurious factors, which need to be
705 ;;; filtered out). Reference: Kenneth Weber, "The accelerated integer
706 ;;; GCD algorithm", ACM Transactions on Mathematical Software, volume
707 ;;; 21, number 1, March 1995, epp. 111-122.
708 (defun bignum-gcd (u0 v0)
709   (declare (type bignum-type u0 v0))
710   (let* ((u1 (if (%bignum-0-or-plusp u0 (%bignum-length u0))
711                  u0
712                  (negate-bignum u0 nil)))
713          (v1 (if (%bignum-0-or-plusp v0 (%bignum-length v0))
714                  v0
715                  (negate-bignum v0 nil))))
716     (if (zerop v1)
717         (return-from bignum-gcd u1))
718     (when (> u1 v1)
719       (rotatef u1 v1))
720     (let ((n (mod v1 u1)))
721       (setf v1 (if (fixnump n)
722                    (make-small-bignum n)
723                    n)))
724     (if (and (= 1 (%bignum-length v1))
725              (zerop (%bignum-ref v1 0)))
726         (return-from bignum-gcd (%normalize-bignum u1
727                                                    (%bignum-length u1))))
728     (let* ((buffer-len (+ 2 (%bignum-length u1)))
729            (u (%allocate-bignum buffer-len))
730            (u-len (%bignum-length u1))
731            (v (%allocate-bignum buffer-len))
732            (v-len (%bignum-length v1))
733            (tmp1 (%allocate-bignum buffer-len))
734            (tmp1-len 0)
735            (tmp2 (%allocate-bignum buffer-len))
736            (tmp2-len 0)
737            (factors-of-two
738             (bignum-factors-of-two u1 (%bignum-length u1)
739                                    v1 (%bignum-length v1))))
740       (declare (type (or null bignum-index)
741                      buffer-len u-len v-len tmp1-len tmp2-len))
742       (bignum-replace u u1)
743       (bignum-replace v v1)
744       (setf u-len
745             (make-gcd-bignum-odd u
746                                  (bignum-buffer-ashift-right u u-len
747                                                              factors-of-two)))
748       (setf v-len
749             (make-gcd-bignum-odd v
750                                  (bignum-buffer-ashift-right v v-len
751                                                              factors-of-two)))
752       (loop until (or (< u-len *accelerated-gcd-cutoff*)
753                       (not v-len)
754                       (zerop v-len)
755                       (and (= 1 v-len)
756                            (zerop (%bignum-ref v 0))))
757         do
758         (gcd-assert (= buffer-len (%bignum-length u)
759                        (%bignum-length v)
760                        (%bignum-length tmp1)
761                        (%bignum-length tmp2)))
762         (if (> (bignum-buffer-integer-length u u-len)
763                (+ #.(truncate sb!vm:n-word-bits 4)
764                   (bignum-buffer-integer-length v v-len)))
765             (setf u-len (dmod u u-len
766                               v v-len
767                               tmp1))
768             (multiple-value-bind (n d) (reduced-ratio-mod u v)
769               (setf tmp1-len
770                     (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
771                                                                    n tmp1))
772               (setf tmp2-len
773                     (multiply-bignum-buffer-and-smallnum-to-buffer u u-len
774                                                                    d tmp2))
775               (gcd-assert (= (copy-bignum tmp2 tmp2-len)
776                              (* (copy-bignum u u-len) d)))
777               (gcd-assert (= (copy-bignum tmp1 tmp1-len)
778                              (* (copy-bignum v v-len) n)))
779               (setf u-len
780                     (subtract-bignum-buffers-with-len tmp1 tmp1-len
781                                                       tmp2 tmp2-len
782                                                       u
783                                                       (1+ (max tmp1-len
784                                                                tmp2-len))))
785               (gcd-assert (or (zerop (- (copy-bignum tmp1 tmp1-len)
786                                         (copy-bignum tmp2 tmp2-len)))
787                               (= (copy-bignum u u-len)
788                                  (- (copy-bignum tmp1 tmp1-len)
789                                     (copy-bignum tmp2 tmp2-len)))))
790               (bignum-abs-buffer u u-len)
791               (gcd-assert (zerop (modularly u)))))
792         (setf u-len (make-gcd-bignum-odd u u-len))
793         (rotatef u v)
794         (rotatef u-len v-len))
795       (setf u (copy-bignum u u-len))
796       (let ((n (bignum-mod-gcd v1 u)))
797         (ash (bignum-mod-gcd u1 (if (fixnump n)
798                                     (make-small-bignum n)
799                                     n))
800              factors-of-two)))))
801
802 (defun bignum-mod-gcd (a b)
803   (declare (type bignum-type a b))
804   (when (< a b)
805     (rotatef a b))
806   ;; While the length difference of A and B is sufficiently large,
807   ;; reduce using MOD (slowish, but it should equalize the sizes of
808   ;; A and B pretty quickly). After that, use the binary GCD
809   ;; algorithm to handle the rest.
810   (loop until (and (= (%bignum-length b) 1) (zerop (%bignum-ref b 0))) do
811         (when (<= (%bignum-length a) (1+ (%bignum-length b)))
812           (return-from bignum-mod-gcd (bignum-binary-gcd a b)))
813         (let ((rem (mod a b)))
814           (if (fixnump rem)
815               (setf a (make-small-bignum rem))
816               (setf a rem))
817           (rotatef a b)))
818   (if (= (%bignum-length a) 1)
819       (%normalize-bignum a 1)
820       a))
821
822 (defun bignum-binary-gcd (a b)
823   (declare (type bignum-type a b))
824   (let* ((len-a (%bignum-length a))
825          (len-b (%bignum-length b)))
826     (declare (type bignum-index len-a len-b))
827     (with-bignum-buffers ((a-buffer len-a a)
828                           (b-buffer len-b b)
829                           (res-buffer (max len-a len-b)))
830       (let* ((factors-of-two
831               (bignum-factors-of-two a-buffer len-a
832                                      b-buffer len-b))
833              (len-a (make-gcd-bignum-odd
834                      a-buffer
835                      (bignum-buffer-ashift-right a-buffer len-a
836                                                  factors-of-two)))
837              (len-b (make-gcd-bignum-odd
838                      b-buffer
839                      (bignum-buffer-ashift-right b-buffer len-b
840                                                  factors-of-two))))
841         (declare (type bignum-index len-a len-b))
842         (let ((x a-buffer)
843               (len-x len-a)
844               (y b-buffer)
845               (len-y len-b)
846               (z res-buffer))
847           (loop
848             (multiple-value-bind (u v len-v r len-r)
849                 (bignum-gcd-order-and-subtract x len-x y len-y z)
850               (declare (type bignum-index len-v len-r))
851               (when (and (= len-r 1) (zerop (%bignum-ref r 0)))
852                 (if (zerop factors-of-two)
853                     (let ((ret (%allocate-bignum len-v)))
854                       (dotimes (i len-v)
855                         (setf (%bignum-ref ret i) (%bignum-ref v i)))
856                       (return (%normalize-bignum ret len-v)))
857                     (return (bignum-ashift-left v factors-of-two len-v))))
858               (setf x v  len-x len-v)
859               (setf y r  len-y (make-gcd-bignum-odd r len-r))
860               (setf z u))))))))
861
862 (defun bignum-gcd-order-and-subtract (a len-a b len-b res)
863   (declare (type bignum-index len-a len-b) (type bignum-type a b))
864   (cond ((= len-a len-b)
865          (do ((i (1- len-a) (1- i)))
866              ((= i -1)
867               (setf (%bignum-ref res 0) 0)
868               (values a b len-b res 1))
869            (let ((a-digit (%bignum-ref a i))
870                  (b-digit (%bignum-ref b i)))
871              (cond ((%digit-compare a-digit b-digit))
872                    ((%digit-greater a-digit b-digit)
873                     (return
874                      (values a b len-b res
875                              (subtract-bignum-buffers a len-a b len-b
876                                                       res))))
877                    (t
878                     (return
879                      (values b a len-a res
880                              (subtract-bignum-buffers b len-b
881                                                       a len-a
882                                                       res))))))))
883         ((> len-a len-b)
884          (values a b len-b res
885                  (subtract-bignum-buffers a len-a b len-b res)))
886         (t
887          (values b a len-a res
888                  (subtract-bignum-buffers b len-b a len-a res)))))
889
890 (defun make-gcd-bignum-odd (a len-a)
891   (declare (type bignum-type a) (type bignum-index len-a))
892   (dotimes (index len-a)
893     (declare (type bignum-index index))
894     (do ((digit (%bignum-ref a index) (%ashr digit 1))
895          (increment 0 (1+ increment)))
896         ((zerop digit))
897       (declare (type (mod #.sb!vm:n-word-bits) increment))
898       (when (oddp digit)
899         (return-from make-gcd-bignum-odd
900                      (bignum-buffer-ashift-right a len-a
901                                                  (+ (* index digit-size)
902                                                     increment)))))))
903
904 \f
905 ;;;; negation
906
907 (eval-when (:compile-toplevel :execute)
908
909 ;;; This negates bignum-len digits of bignum, storing the resulting digits into
910 ;;; result (possibly EQ to bignum) and returning whatever end-carry there is.
911 (sb!xc:defmacro bignum-negate-loop (bignum
912                                     bignum-len
913                                     &optional (result nil resultp))
914   (let ((carry (gensym))
915         (end (gensym))
916         (value (gensym))
917         (last (gensym)))
918     `(let* (,@(if (not resultp) `(,last))
919             (,carry
920              (multiple-value-bind (,value ,carry)
921                  (%add-with-carry (%lognot (%bignum-ref ,bignum 0)) 1 0)
922                ,(if resultp
923                     `(setf (%bignum-ref ,result 0) ,value)
924                     `(setf ,last ,value))
925                ,carry))
926             (i 1)
927             (,end ,bignum-len))
928        (declare (type bit ,carry)
929                 (type bignum-index i ,end))
930        (loop
931          (when (= i ,end) (return))
932          (multiple-value-bind (,value temp)
933              (%add-with-carry (%lognot (%bignum-ref ,bignum i)) 0 ,carry)
934            ,(if resultp
935                 `(setf (%bignum-ref ,result i) ,value)
936                 `(setf ,last ,value))
937            (setf ,carry temp))
938          (incf i))
939        ,(if resultp carry `(values ,carry ,last)))))
940
941 ) ; EVAL-WHEN
942
943 ;;; Fully-normalize is an internal optional. It cause this to always return
944 ;;; a bignum, without any extraneous digits, and it never returns a fixnum.
945 (defun negate-bignum (x &optional (fully-normalize t))
946   (declare (type bignum-type x))
947   (let* ((len-x (%bignum-length x))
948          (len-res (1+ len-x))
949          (res (%allocate-bignum len-res)))
950     (declare (type bignum-index len-x len-res)) ;Test len-res for range?
951     (let ((carry (bignum-negate-loop x len-x res)))
952       (setf (%bignum-ref res len-x)
953             (%add-with-carry (%lognot (%sign-digit x len-x)) 0 carry)))
954     (if fully-normalize
955         (%normalize-bignum res len-res)
956         (%mostly-normalize-bignum res len-res))))
957
958 ;;; This assumes bignum is positive; that is, the result of negating it will
959 ;;; stay in the provided allocated bignum.
960 (defun negate-bignum-buffer-in-place (bignum bignum-len)
961   (bignum-negate-loop bignum bignum-len bignum)
962   bignum)
963
964 (defun negate-bignum-in-place (bignum)
965   (declare (inline negate-bignum-buffer-in-place))
966   (negate-bignum-buffer-in-place bignum (%bignum-length bignum)))
967
968 (defun bignum-abs-buffer (bignum len)
969   (unless (%bignum-0-or-plusp bignum len)
970     (negate-bignum-buffer-in-place bignum len)))
971 \f
972 ;;;; shifting
973
974 (eval-when (:compile-toplevel :execute)
975
976 ;;; This macro is used by BIGNUM-ASHIFT-RIGHT, BIGNUM-BUFFER-ASHIFT-RIGHT, and
977 ;;; BIGNUM-LDB-BIGNUM-RES. They supply a termination form that references
978 ;;; locals established by this form. Source is the source bignum. Start-digit
979 ;;; is the first digit in source from which we pull bits. Start-pos is the
980 ;;; first bit we want. Res-len-form is the form that computes the length of
981 ;;; the resulting bignum. Termination is a DO termination form with a test and
982 ;;; body. When result is supplied, it is the variable to which this binds a
983 ;;; newly allocated bignum.
984 ;;;
985 ;;; Given start-pos, 1-31 inclusively, of shift, we form the j'th resulting
986 ;;; digit from high bits of the i'th source digit and the start-pos number of
987 ;;; bits from the i+1'th source digit.
988 (sb!xc:defmacro shift-right-unaligned (source
989                                        start-digit
990                                        start-pos
991                                        res-len-form
992                                        termination
993                                        &optional result)
994   `(let* ((high-bits-in-first-digit (- digit-size ,start-pos))
995           (res-len ,res-len-form)
996           (res-len-1 (1- res-len))
997           ,@(if result `((,result (%allocate-bignum res-len)))))
998      (declare (type bignum-index res-len res-len-1))
999      (do ((i ,start-digit (1+ i))
1000           (j 0 (1+ j)))
1001          ,termination
1002        (declare (type bignum-index i j))
1003        (setf (%bignum-ref ,(if result result source) j)
1004              (%logior (%digit-logical-shift-right (%bignum-ref ,source i)
1005                                                   ,start-pos)
1006                       (%ashl (%bignum-ref ,source (1+ i))
1007                              high-bits-in-first-digit))))))
1008
1009 ) ; EVAL-WHEN
1010
1011 ;;; First compute the number of whole digits to shift, shifting them by
1012 ;;; skipping them when we start to pick up bits, and the number of bits to
1013 ;;; shift the remaining digits into place. If the number of digits is greater
1014 ;;; than the length of the bignum, then the result is either 0 or -1. If we
1015 ;;; shift on a digit boundary (that is, n-bits is zero), then we just copy
1016 ;;; digits. The last branch handles the general case which uses a macro that a
1017 ;;; couple other routines use. The fifth argument to the macro references
1018 ;;; locals established by the macro.
1019 (defun bignum-ashift-right (bignum count)
1020   (declare (type bignum-type bignum)
1021            (type unsigned-byte count))
1022   (let ((bignum-len (%bignum-length bignum)))
1023     (declare (type bignum-index bignum-len))
1024     (cond ((fixnump count)
1025            (multiple-value-bind (digits n-bits) (truncate count digit-size)
1026              (declare (type bignum-index digits))
1027              (cond
1028               ((>= digits bignum-len)
1029                (if (%bignum-0-or-plusp bignum bignum-len) 0 -1))
1030               ((zerop n-bits)
1031                (bignum-ashift-right-digits bignum digits))
1032               (t
1033                (shift-right-unaligned bignum digits n-bits (- bignum-len digits)
1034                                       ((= j res-len-1)
1035                                        (setf (%bignum-ref res j)
1036                                              (%ashr (%bignum-ref bignum i) n-bits))
1037                                        (%normalize-bignum res res-len))
1038                                       res)))))
1039           ((> count bignum-len)
1040            (if (%bignum-0-or-plusp bignum bignum-len) 0 -1))
1041            ;; Since a FIXNUM should be big enough to address anything in
1042            ;; memory, including arrays of bits, and since arrays of bits
1043            ;; take up about the same space as corresponding fixnums, there
1044            ;; should be no way that we fall through to this case: any shift
1045            ;; right by a bignum should give zero. But let's check anyway:
1046           (t (error "bignum overflow: can't shift right by ~S" count)))))
1047
1048 (defun bignum-ashift-right-digits (bignum digits)
1049   (declare (type bignum-type bignum)
1050            (type bignum-index digits))
1051   (let* ((res-len (- (%bignum-length bignum) digits))
1052          (res (%allocate-bignum res-len)))
1053     (declare (type bignum-index res-len)
1054              (type bignum-type res))
1055     (bignum-replace res bignum :start2 digits)
1056     (%normalize-bignum res res-len)))
1057
1058 ;;; GCD uses this for an in-place shifting operation. This is different enough
1059 ;;; from BIGNUM-ASHIFT-RIGHT that it isn't worth folding the bodies into a
1060 ;;; macro, but they share the basic algorithm. This routine foregoes a first
1061 ;;; test for digits being greater than or equal to bignum-len since that will
1062 ;;; never happen for its uses in GCD. We did fold the last branch into a macro
1063 ;;; since it was duplicated a few times, and the fifth argument to it
1064 ;;; references locals established by the macro.
1065 (defun bignum-buffer-ashift-right (bignum bignum-len x)
1066   (declare (type bignum-index bignum-len) (fixnum x))
1067   (multiple-value-bind (digits n-bits) (truncate x digit-size)
1068     (declare (type bignum-index digits))
1069     (cond
1070      ((zerop n-bits)
1071       (let ((new-end (- bignum-len digits)))
1072         (bignum-replace bignum bignum :end1 new-end :start2 digits
1073                         :end2 bignum-len)
1074         (%normalize-bignum-buffer bignum new-end)))
1075      (t
1076       (shift-right-unaligned bignum digits n-bits (- bignum-len digits)
1077                              ((= j res-len-1)
1078                               (setf (%bignum-ref bignum j)
1079                                     (%ashr (%bignum-ref bignum i) n-bits))
1080                               (%normalize-bignum-buffer bignum res-len)))))))
1081
1082 ;;; This handles shifting a bignum buffer to provide fresh bignum data for some
1083 ;;; internal routines. We know bignum is safe when called with bignum-len.
1084 ;;; First we compute the number of whole digits to shift, shifting them
1085 ;;; starting to store farther along the result bignum. If we shift on a digit
1086 ;;; boundary (that is, n-bits is zero), then we just copy digits. The last
1087 ;;; branch handles the general case.
1088 (defun bignum-ashift-left (bignum x &optional bignum-len)
1089   (declare (type bignum-type bignum)
1090            (type unsigned-byte x)
1091            (type (or null bignum-index) bignum-len))
1092   (if (fixnump x)
1093     (multiple-value-bind (digits n-bits) (truncate x digit-size)
1094       (let* ((bignum-len (or bignum-len (%bignum-length bignum)))
1095              (res-len (+ digits bignum-len 1)))
1096         (when (> res-len maximum-bignum-length)
1097           (error "can't represent result of left shift"))
1098         (if (zerop n-bits)
1099           (bignum-ashift-left-digits bignum bignum-len digits)
1100           (bignum-ashift-left-unaligned bignum digits n-bits res-len))))
1101     ;; Left shift by a number too big to be represented as a fixnum
1102     ;; would exceed our memory capacity, since a fixnum is big enough
1103     ;; to index any array, including a bit array.
1104     (error "can't represent result of left shift")))
1105
1106 (defun bignum-ashift-left-digits (bignum bignum-len digits)
1107   (declare (type bignum-index bignum-len digits))
1108   (let* ((res-len (+ bignum-len digits))
1109          (res (%allocate-bignum res-len)))
1110     (declare (type bignum-index res-len))
1111     (bignum-replace res bignum :start1 digits :end1 res-len :end2 bignum-len
1112                     :from-end t)
1113     res))
1114
1115 ;;; BIGNUM-TRUNCATE uses this to store into a bignum buffer by supplying res.
1116 ;;; When res comes in non-nil, then this foregoes allocating a result, and it
1117 ;;; normalizes the buffer instead of the would-be allocated result.
1118 ;;;
1119 ;;; We start storing into one digit higher than digits, storing a whole result
1120 ;;; digit from parts of two contiguous digits from bignum. When the loop
1121 ;;; finishes, we store the remaining bits from bignum's first digit in the
1122 ;;; first non-zero result digit, digits. We also grab some left over high
1123 ;;; bits from the last digit of bignum.
1124 (defun bignum-ashift-left-unaligned (bignum digits n-bits res-len
1125                                      &optional (res nil resp))
1126   (declare (type bignum-index digits res-len)
1127            (type (mod #.digit-size) n-bits))
1128   (let* ((remaining-bits (- digit-size n-bits))
1129          (res-len-1 (1- res-len))
1130          (res (or res (%allocate-bignum res-len))))
1131     (declare (type bignum-index res-len res-len-1))
1132     (do ((i 0 (1+ i))
1133          (j (1+ digits) (1+ j)))
1134         ((= j res-len-1)
1135          (setf (%bignum-ref res digits)
1136                (%ashl (%bignum-ref bignum 0) n-bits))
1137          (setf (%bignum-ref res j)
1138                (%ashr (%bignum-ref bignum i) remaining-bits))
1139          (if resp
1140              (%normalize-bignum-buffer res res-len)
1141              (%normalize-bignum res res-len)))
1142       (declare (type bignum-index i j))
1143       (setf (%bignum-ref res j)
1144             (%logior (%digit-logical-shift-right (%bignum-ref bignum i)
1145                                                  remaining-bits)
1146                      (%ashl (%bignum-ref bignum (1+ i)) n-bits))))))
1147 \f
1148 ;;;; relational operators
1149
1150 ;;; Return T iff bignum is positive.
1151 (defun bignum-plus-p (bignum)
1152   (declare (type bignum-type bignum))
1153   (%bignum-0-or-plusp bignum (%bignum-length bignum)))
1154
1155 ;;; This compares two bignums returning -1, 0, or 1, depending on
1156 ;;; whether a is less than, equal to, or greater than b.
1157 (declaim (ftype (function (bignum bignum) (integer -1 1)) bignum-compare))
1158 (defun bignum-compare (a b)
1159   (declare (type bignum-type a b))
1160   (let* ((len-a (%bignum-length a))
1161          (len-b (%bignum-length b))
1162          (a-plusp (%bignum-0-or-plusp a len-a))
1163          (b-plusp (%bignum-0-or-plusp b len-b)))
1164     (declare (type bignum-index len-a len-b))
1165     (cond ((not (eq a-plusp b-plusp))
1166            (if a-plusp 1 -1))
1167           ((= len-a len-b)
1168            (do ((i (1- len-a) (1- i)))
1169                (())
1170              (declare (type bignum-index i))
1171              (let ((a-digit (%bignum-ref a i))
1172                    (b-digit (%bignum-ref b i)))
1173                (declare (type bignum-element-type a-digit b-digit))
1174                (when (%digit-greater a-digit b-digit)
1175                  (return 1))
1176                (when (%digit-greater b-digit a-digit)
1177                  (return -1)))
1178              (when (zerop i) (return 0))))
1179           ((> len-a len-b)
1180            (if a-plusp 1 -1))
1181           (t (if a-plusp -1 1)))))
1182 \f
1183 ;;;; float conversion
1184
1185 ;;; Make a single or double float with the specified significand,
1186 ;;; exponent and sign.
1187 (defun single-float-from-bits (bits exp plusp)
1188   (declare (fixnum exp))
1189   (declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
1190   (let ((res (dpb exp
1191                   sb!vm:single-float-exponent-byte
1192                   (logandc2 (logand #xffffffff
1193                                     (%bignum-ref bits 1))
1194                             sb!vm:single-float-hidden-bit))))
1195     (make-single-float
1196      (if plusp
1197          res
1198          (logior res (ash -1 sb!vm:float-sign-shift))))))
1199 (defun double-float-from-bits (bits exp plusp)
1200   (declare (fixnum exp))
1201   (declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
1202   (let ((hi (dpb exp
1203                  sb!vm:double-float-exponent-byte
1204                  (logandc2 (ecase sb!vm::n-word-bits
1205                              (32 (%bignum-ref bits 2))
1206                              (64 (ash (%bignum-ref bits 1) -32)))
1207                            sb!vm:double-float-hidden-bit)))
1208         (lo (logand #xffffffff (%bignum-ref bits 1))))
1209     (make-double-float (if plusp
1210                            hi
1211                            (logior hi (ash -1 sb!vm:float-sign-shift)))
1212                        lo)))
1213 #!+(and long-float x86)
1214 (defun long-float-from-bits (bits exp plusp)
1215   (declare (fixnum exp))
1216   (declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
1217   (make-long-float
1218    (if plusp
1219        exp
1220        (logior exp (ash 1 15)))
1221    (%bignum-ref bits 2)
1222    (%bignum-ref bits 1)))
1223
1224 ;;; Convert Bignum to a float in the specified Format, rounding to the best
1225 ;;; approximation.
1226 (defun bignum-to-float (bignum format)
1227   (let* ((plusp (bignum-plus-p bignum))
1228          (x (if plusp bignum (negate-bignum bignum)))
1229          (len (bignum-integer-length x))
1230          (digits (float-format-digits format))
1231          (keep (+ digits digit-size))
1232          (shift (- keep len))
1233          (shifted (if (minusp shift)
1234                       (bignum-ashift-right x (- shift))
1235                       (bignum-ashift-left x shift)))
1236          (low (%bignum-ref shifted 0))
1237          (round-bit (ash 1 (1- digit-size))))
1238     (declare (type bignum-index len digits keep) (fixnum shift))
1239     (labels ((round-up ()
1240                (let ((rounded (add-bignums shifted round-bit)))
1241                  (if (> (integer-length rounded) keep)
1242                      (float-from-bits (bignum-ashift-right rounded 1)
1243                                       (1+ len))
1244                      (float-from-bits rounded len))))
1245              (float-from-bits (bits len)
1246                (declare (type bignum-index len))
1247                (ecase format
1248                  (single-float
1249                   (single-float-from-bits
1250                    bits
1251                    (check-exponent len sb!vm:single-float-bias
1252                                    sb!vm:single-float-normal-exponent-max)
1253                    plusp))
1254                  (double-float
1255                   (double-float-from-bits
1256                    bits
1257                    (check-exponent len sb!vm:double-float-bias
1258                                    sb!vm:double-float-normal-exponent-max)
1259                    plusp))
1260                  #!+long-float
1261                  (long-float
1262                   (long-float-from-bits
1263                    bits
1264                    (check-exponent len sb!vm:long-float-bias
1265                                    sb!vm:long-float-normal-exponent-max)
1266                    plusp))))
1267              (check-exponent (exp bias max)
1268                (declare (type bignum-index len))
1269                (let ((exp (+ exp bias)))
1270                  (when (> exp max)
1271                    ;; Why a SIMPLE-TYPE-ERROR? Well, this is mainly
1272                    ;; called by COERCE, which requires an error of
1273                    ;; TYPE-ERROR if the conversion can't happen
1274                    ;; (except in certain circumstances when we are
1275                    ;; coercing to a FUNCTION) -- CSR, 2002-09-18
1276                    (error 'simple-type-error
1277                           :format-control "Too large to be represented as a ~S:~%  ~S"
1278                           :format-arguments (list format x)
1279                           :expected-type format
1280                           :datum x))
1281                  exp)))
1282
1283     (cond
1284      ;; Round down if round bit is 0.
1285      ((not (logtest round-bit low))
1286       (float-from-bits shifted len))
1287      ;; If only round bit is set, then round to even.
1288      ((and (= low round-bit)
1289            (dotimes (i (- (%bignum-length x) (ceiling keep digit-size))
1290                        t)
1291              (unless (zerop (%bignum-ref x i)) (return nil))))
1292       (let ((next (%bignum-ref shifted 1)))
1293         (if (oddp next)
1294             (round-up)
1295             (float-from-bits shifted len))))
1296      ;; Otherwise, round up.
1297      (t
1298       (round-up))))))
1299 \f
1300 ;;;; integer length and logbitp/logcount
1301
1302 (defun bignum-buffer-integer-length (bignum len)
1303   (declare (type bignum-type bignum))
1304   (let* ((len-1 (1- len))
1305          (digit (%bignum-ref bignum len-1)))
1306     (declare (type bignum-index len len-1)
1307              (type bignum-element-type digit))
1308     (+ (integer-length (%fixnum-digit-with-correct-sign digit))
1309        (* len-1 digit-size))))
1310
1311 (defun bignum-integer-length (bignum)
1312   (declare (type bignum-type bignum))
1313   (bignum-buffer-integer-length bignum (%bignum-length bignum)))
1314
1315 (defun bignum-logbitp (index bignum)
1316   (declare (type bignum-type bignum))
1317   (let ((len (%bignum-length bignum)))
1318     (declare (type bignum-index len))
1319     (multiple-value-bind (word-index bit-index)
1320         (floor index digit-size)
1321       (if (>= word-index len)
1322           (not (bignum-plus-p bignum))
1323           (logbitp bit-index (%bignum-ref bignum word-index))))))
1324
1325 (defun bignum-logcount (bignum)
1326   (declare (type bignum-type bignum))
1327   (let ((length (%bignum-length bignum))
1328         (result 0))
1329     (declare (type bignum-index length)
1330              (fixnum result))
1331     (do ((index 0 (1+ index)))
1332         ((= index length)
1333          (if (%bignum-0-or-plusp bignum length)
1334              result
1335              (- (* length digit-size) result)))
1336       (let ((digit (%bignum-ref bignum index)))
1337         (declare (type bignum-element-type digit))
1338         (incf result (logcount digit))))))
1339 \f
1340 ;;;; logical operations
1341
1342 ;;;; NOT
1343
1344 (defun bignum-logical-not (a)
1345   (declare (type bignum-type a))
1346   (let* ((len (%bignum-length a))
1347          (res (%allocate-bignum len)))
1348     (declare (type bignum-index len))
1349     (dotimes (i len res)
1350       (declare (type bignum-index i))
1351       (setf (%bignum-ref res i) (%lognot (%bignum-ref a i))))))
1352
1353 ;;;; AND
1354
1355 (defun bignum-logical-and (a b)
1356   (declare (type bignum-type a b))
1357   (let* ((len-a (%bignum-length a))
1358          (len-b (%bignum-length b))
1359          (a-plusp (%bignum-0-or-plusp a len-a))
1360          (b-plusp (%bignum-0-or-plusp b len-b)))
1361     (declare (type bignum-index len-a len-b))
1362     (cond
1363      ((< len-a len-b)
1364       (if a-plusp
1365           (logand-shorter-positive a len-a b (%allocate-bignum len-a))
1366           (logand-shorter-negative a len-a b len-b (%allocate-bignum len-b))))
1367      ((< len-b len-a)
1368       (if b-plusp
1369           (logand-shorter-positive b len-b a (%allocate-bignum len-b))
1370           (logand-shorter-negative b len-b a len-a (%allocate-bignum len-a))))
1371      (t (logand-shorter-positive a len-a b (%allocate-bignum len-a))))))
1372
1373 ;;; This takes a shorter bignum, a and len-a, that is positive. Because this
1374 ;;; is AND, we don't care about any bits longer than a's since its infinite 0
1375 ;;; sign bits will mask the other bits out of b. The result is len-a big.
1376 (defun logand-shorter-positive (a len-a b res)
1377   (declare (type bignum-type a b res)
1378            (type bignum-index len-a))
1379   (dotimes (i len-a)
1380     (declare (type bignum-index i))
1381     (setf (%bignum-ref res i)
1382           (%logand (%bignum-ref a i) (%bignum-ref b i))))
1383   (%normalize-bignum res len-a))
1384
1385 ;;; This takes a shorter bignum, a and len-a, that is negative. Because this
1386 ;;; is AND, we just copy any bits longer than a's since its infinite 1 sign
1387 ;;; bits will include any bits from b. The result is len-b big.
1388 (defun logand-shorter-negative (a len-a b len-b res)
1389   (declare (type bignum-type a b res)
1390            (type bignum-index len-a len-b))
1391   (dotimes (i len-a)
1392     (declare (type bignum-index i))
1393     (setf (%bignum-ref res i)
1394           (%logand (%bignum-ref a i) (%bignum-ref b i))))
1395   (do ((i len-a (1+ i)))
1396       ((= i len-b))
1397     (declare (type bignum-index i))
1398     (setf (%bignum-ref res i) (%bignum-ref b i)))
1399   (%normalize-bignum res len-b))
1400
1401 ;;;; IOR
1402
1403 (defun bignum-logical-ior (a b)
1404   (declare (type bignum-type a b))
1405   (let* ((len-a (%bignum-length a))
1406          (len-b (%bignum-length b))
1407          (a-plusp (%bignum-0-or-plusp a len-a))
1408          (b-plusp (%bignum-0-or-plusp b len-b)))
1409     (declare (type bignum-index len-a len-b))
1410     (cond
1411      ((< len-a len-b)
1412       (if a-plusp
1413           (logior-shorter-positive a len-a b len-b (%allocate-bignum len-b))
1414           (logior-shorter-negative a len-a b len-b (%allocate-bignum len-b))))
1415      ((< len-b len-a)
1416       (if b-plusp
1417           (logior-shorter-positive b len-b a len-a (%allocate-bignum len-a))
1418           (logior-shorter-negative b len-b a len-a (%allocate-bignum len-a))))
1419      (t (logior-shorter-positive a len-a b len-b (%allocate-bignum len-a))))))
1420
1421 ;;; This takes a shorter bignum, a and len-a, that is positive. Because this
1422 ;;; is IOR, we don't care about any bits longer than a's since its infinite
1423 ;;; 0 sign bits will mask the other bits out of b out to len-b. The result
1424 ;;; is len-b long.
1425 (defun logior-shorter-positive (a len-a b len-b res)
1426   (declare (type bignum-type a b res)
1427            (type bignum-index len-a len-b))
1428   (dotimes (i len-a)
1429     (declare (type bignum-index i))
1430     (setf (%bignum-ref res i)
1431           (%logior (%bignum-ref a i) (%bignum-ref b i))))
1432   (do ((i len-a (1+ i)))
1433       ((= i len-b))
1434     (declare (type bignum-index i))
1435     (setf (%bignum-ref res i) (%bignum-ref b i)))
1436   (%normalize-bignum res len-b))
1437
1438 ;;; This takes a shorter bignum, a and len-a, that is negative. Because this
1439 ;;; is IOR, we just copy any bits longer than a's since its infinite 1 sign
1440 ;;; bits will include any bits from b. The result is len-b long.
1441 (defun logior-shorter-negative (a len-a b len-b res)
1442   (declare (type bignum-type a b res)
1443            (type bignum-index len-a len-b))
1444   (dotimes (i len-a)
1445     (declare (type bignum-index i))
1446     (setf (%bignum-ref res i)
1447           (%logior (%bignum-ref a i) (%bignum-ref b i))))
1448   (do ((i len-a (1+ i))
1449        (sign (%sign-digit a len-a)))
1450       ((= i len-b))
1451     (declare (type bignum-index i))
1452     (setf (%bignum-ref res i) sign))
1453   (%normalize-bignum res len-b))
1454
1455 ;;;; XOR
1456
1457 (defun bignum-logical-xor (a b)
1458   (declare (type bignum-type a b))
1459   (let ((len-a (%bignum-length a))
1460         (len-b (%bignum-length b)))
1461     (declare (type bignum-index len-a len-b))
1462     (if (< len-a len-b)
1463         (bignum-logical-xor-aux a len-a b len-b (%allocate-bignum len-b))
1464         (bignum-logical-xor-aux b len-b a len-a (%allocate-bignum len-a)))))
1465
1466 ;;; This takes the shorter of two bignums in a and len-a. Res is len-b
1467 ;;; long. Do the XOR.
1468 (defun bignum-logical-xor-aux (a len-a b len-b res)
1469   (declare (type bignum-type a b res)
1470            (type bignum-index len-a len-b))
1471   (dotimes (i len-a)
1472     (declare (type bignum-index i))
1473     (setf (%bignum-ref res i)
1474           (%logxor (%bignum-ref a i) (%bignum-ref b i))))
1475   (do ((i len-a (1+ i))
1476        (sign (%sign-digit a len-a)))
1477       ((= i len-b))
1478     (declare (type bignum-index i))
1479     (setf (%bignum-ref res i) (%logxor sign (%bignum-ref b i))))
1480   (%normalize-bignum res len-b))
1481 \f
1482 ;;;; There used to be a bunch of code to implement "efficient" versions of LDB
1483 ;;;; and DPB here.  But it apparently was never used, so it's been deleted.
1484 ;;;;   --njf, 2007-02-04
1485 \f
1486 ;;;; TRUNCATE
1487
1488 ;;; This is the original sketch of the algorithm from which I implemented this
1489 ;;; TRUNCATE, assuming both operands are bignums. I should modify this to work
1490 ;;; with the documentation on my functions, as a general introduction. I've
1491 ;;; left this here just in case someone needs it in the future. Don't look at
1492 ;;; this unless reading the functions' comments leaves you at a loss. Remember
1493 ;;; this comes from Knuth, so the book might give you the right general
1494 ;;; overview.
1495 ;;;
1496 ;;; (truncate x y):
1497 ;;;
1498 ;;; If X's magnitude is less than Y's, then result is 0 with remainder X.
1499 ;;;
1500 ;;; Make x and y positive, copying x if it is already positive.
1501 ;;;
1502 ;;; Shift y left until there's a 1 in the 30'th bit (most significant, non-sign
1503 ;;;       digit)
1504 ;;;    Just do most sig digit to determine how much to shift whole number.
1505 ;;; Shift x this much too.
1506 ;;; Remember this initial shift count.
1507 ;;;
1508 ;;; Allocate q to be len-x minus len-y quantity plus 1.
1509 ;;;
1510 ;;; i = last digit of x.
1511 ;;; k = last digit of q.
1512 ;;;
1513 ;;; LOOP
1514 ;;;
1515 ;;; j = last digit of y.
1516 ;;;
1517 ;;; compute guess.
1518 ;;; if x[i] = y[j] then g = (1- (ash 1 digit-size))
1519 ;;; else g = x[i]x[i-1]/y[j].
1520 ;;;
1521 ;;; check guess.
1522 ;;; %UNSIGNED-MULTIPLY returns b and c defined below.
1523 ;;;    a = x[i-1] - (logand (* g y[j]) #xFFFFFFFF).
1524 ;;;       Use %UNSIGNED-MULTIPLY taking low-order result.
1525 ;;;    b = (logand (ash (* g y[j-1]) (- digit-size)) (1- (ash 1 digit-size))).
1526 ;;;    c = (logand (* g y[j-1]) (1- (ash 1 digit-size))).
1527 ;;; if a < b, okay.
1528 ;;; if a > b, guess is too high
1529 ;;;    g = g - 1; go back to "check guess".
1530 ;;; if a = b and c > x[i-2], guess is too high
1531 ;;;    g = g - 1; go back to "check guess".
1532 ;;; GUESS IS 32-BIT NUMBER, SO USE THING TO KEEP IN SPECIAL REGISTER
1533 ;;; SAME FOR A, B, AND C.
1534 ;;;
1535 ;;; Subtract g * y from x[i - len-y+1]..x[i]. See paper for doing this in step.
1536 ;;; If x[i] < 0, guess is screwed up.
1537 ;;;    negative g, then add 1
1538 ;;;    zero or positive g, then subtract 1
1539 ;;; AND add y back into x[len-y+1..i].
1540 ;;;
1541 ;;; q[k] = g.
1542 ;;; i = i - 1.
1543 ;;; k = k - 1.
1544 ;;;
1545 ;;; If k>=0, goto LOOP.
1546 ;;;
1547 ;;; Now quotient is good, but remainder is not.
1548 ;;; Shift x right by saved initial left shifting count.
1549 ;;;
1550 ;;; Check quotient and remainder signs.
1551 ;;; x pos y pos --> q pos r pos
1552 ;;; x pos y neg --> q neg r pos
1553 ;;; x neg y pos --> q neg r neg
1554 ;;; x neg y neg --> q pos r neg
1555 ;;;
1556 ;;; Normalize quotient and remainder. Cons result if necessary.
1557
1558
1559 ;;; This used to be split into multiple functions, which shared state
1560 ;;; in special variables *TRUNCATE-X* and *TRUNCATE-Y*. Having so many
1561 ;;; special variable accesses in tight inner loops was having a large
1562 ;;; effect on performance, so the helper functions have now been
1563 ;;; refactored into local functions and the special variables into
1564 ;;; lexicals.  There was also a lot of boxing and unboxing of
1565 ;;; (UNSIGNED-BYTE 32)'s going on, which this refactoring
1566 ;;; eliminated. This improves the performance on some CL-BENCH tests
1567 ;;; by up to 50%, which is probably signigicant enough to justify the
1568 ;;; reduction in readability that was introduced. --JES, 2004-08-07
1569 (defun bignum-truncate (x y)
1570   (declare (type bignum-type x y))
1571   (let (truncate-x truncate-y)
1572     (labels
1573         ;;; Divide X by Y when Y is a single bignum digit. BIGNUM-TRUNCATE
1574         ;;; fixes up the quotient and remainder with respect to sign and
1575         ;;; normalization.
1576         ;;;
1577         ;;; We don't have to worry about shifting Y to make its most
1578         ;;; significant digit sufficiently large for %FLOOR to return
1579         ;;; digit-size quantities for the q-digit and r-digit. If Y is
1580         ;;; a single digit bignum, it is already large enough for
1581         ;;; %FLOOR. That is, it has some bits on pretty high in the
1582         ;;; digit.
1583         ((bignum-truncate-single-digit (x len-x y)
1584            (declare (type bignum-index len-x))
1585            (let ((y (%bignum-ref y 0)))
1586              (declare (type bignum-element-type y))
1587              (if (not (logtest y (1- y)))
1588                  ;; Y is a power of two.
1589                  (if (= y 1)
1590                      ;; SHIFT-RIGHT-UNALIGNED won't do the right thing
1591                      ;; with a shift count of 0, so special case this.
1592                      ;; We could probably get away with (VALUES X 0)
1593                      ;; here, but it's not clear that some of the
1594                      ;; normalization logic further down would avoid
1595                      ;; mutilating X.  Just go ahead and cons, consing's
1596                      ;; cheap.
1597                      (values (copy-bignum x len-x) 0)
1598                      (let ((n-bits (1- (integer-length y))))
1599                        (values
1600                         (shift-right-unaligned x 0 n-bits len-x
1601                                                ((= j res-len-1)
1602                                                 (setf (%bignum-ref res j)
1603                                                       (%ashr (%bignum-ref x i) n-bits))
1604                                                 res)
1605                                                res)
1606                         (logand (%bignum-ref x 0) (1- y)))))
1607                  (do ((i (1- len-x) (1- i))
1608                       (q (%allocate-bignum len-x))
1609                       (r 0))
1610                      ((minusp i)
1611                       (let ((rem (%allocate-bignum 1)))
1612                         (setf (%bignum-ref rem 0) r)
1613                         (values q rem)))
1614                    (declare (type bignum-element-type r))
1615                    (multiple-value-bind (q-digit r-digit)
1616                        (%floor r (%bignum-ref x i) y)
1617                      (declare (type bignum-element-type q-digit r-digit))
1618                      (setf (%bignum-ref q i) q-digit)
1619                      (setf r r-digit))))))
1620         ;;; This returns a guess for the next division step. Y1 is the
1621         ;;; highest y digit, and y2 is the second to highest y
1622         ;;; digit. The x... variables are the three highest x digits
1623         ;;; for the next division step.
1624         ;;;
1625         ;;; From Knuth, our guess is either all ones or x-i and x-i-1
1626         ;;; divided by y1, depending on whether x-i and y1 are the
1627         ;;; same. We test this guess by determining whether guess*y2
1628         ;;; is greater than the three high digits of x minus guess*y1
1629         ;;; shifted left one digit:
1630         ;;;    ------------------------------
1631         ;;;   |    x-i    |   x-i-1  | x-i-2 |
1632         ;;;    ------------------------------
1633         ;;;    ------------------------------
1634         ;;; - | g*y1 high | g*y1 low |   0   |
1635         ;;;    ------------------------------
1636         ;;;             ...               <   guess*y2     ???
1637         ;;; If guess*y2 is greater, then we decrement our guess by one
1638         ;;; and try again.  This returns a guess that is either
1639         ;;; correct or one too large.
1640          (bignum-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
1641            (declare (type bignum-element-type y1 y2 x-i x-i-1 x-i-2))
1642            (let ((guess (if (%digit-compare x-i y1)
1643                             all-ones-digit
1644                             (%floor x-i x-i-1 y1))))
1645              (declare (type bignum-element-type guess))
1646              (loop
1647                  (multiple-value-bind (high-guess*y1 low-guess*y1)
1648                      (%multiply guess y1)
1649                    (declare (type bignum-element-type low-guess*y1
1650                                   high-guess*y1))
1651                    (multiple-value-bind (high-guess*y2 low-guess*y2)
1652                        (%multiply guess y2)
1653                      (declare (type bignum-element-type high-guess*y2
1654                                     low-guess*y2))
1655                      (multiple-value-bind (middle-digit borrow)
1656                          (%subtract-with-borrow x-i-1 low-guess*y1 1)
1657                        (declare (type bignum-element-type middle-digit)
1658                                 (fixnum borrow))
1659                        ;; Supplying borrow of 1 means there was no
1660                        ;; borrow, and we know x-i-2 minus 0 requires
1661                        ;; no borrow.
1662                        (let ((high-digit (%subtract-with-borrow x-i
1663                                                                 high-guess*y1
1664                                                                 borrow)))
1665                          (declare (type bignum-element-type high-digit))
1666                          (if (and (%digit-compare high-digit 0)
1667                                   (or (%digit-greater high-guess*y2
1668                                                       middle-digit)
1669                                       (and (%digit-compare middle-digit
1670                                                            high-guess*y2)
1671                                            (%digit-greater low-guess*y2
1672                                                            x-i-2))))
1673                              (setf guess (%subtract-with-borrow guess 1 1))
1674                              (return guess)))))))))
1675         ;;; Divide TRUNCATE-X by TRUNCATE-Y, returning the quotient
1676         ;;; and destructively modifying TRUNCATE-X so that it holds
1677         ;;; the remainder.
1678         ;;;
1679         ;;; LEN-X and LEN-Y tell us how much of the buffers we care about.
1680         ;;;
1681         ;;; TRUNCATE-X definitely has at least three digits, and it has one
1682         ;;; more than TRUNCATE-Y. This keeps i, i-1, i-2, and low-x-digit
1683         ;;; happy. Thanks to SHIFT-AND-STORE-TRUNCATE-BUFFERS.
1684          (return-quotient-leaving-remainder (len-x len-y)
1685            (declare (type bignum-index len-x len-y))
1686            (let* ((len-q (- len-x len-y))
1687                   ;; Add one for extra sign digit in case high bit is on.
1688                   (q (%allocate-bignum (1+ len-q)))
1689                   (k (1- len-q))
1690                   (y1 (%bignum-ref truncate-y (1- len-y)))
1691                   (y2 (%bignum-ref truncate-y (- len-y 2)))
1692                   (i (1- len-x))
1693                   (i-1 (1- i))
1694                   (i-2 (1- i-1))
1695                   (low-x-digit (- i len-y)))
1696              (declare (type bignum-index len-q k i i-1 i-2 low-x-digit)
1697                       (type bignum-element-type y1 y2))
1698              (loop
1699                  (setf (%bignum-ref q k)
1700                        (try-bignum-truncate-guess
1701                         ;; This modifies TRUNCATE-X. Must access
1702                         ;; elements each pass.
1703                         (bignum-truncate-guess y1 y2
1704                                                (%bignum-ref truncate-x i)
1705                                                (%bignum-ref truncate-x i-1)
1706                                                (%bignum-ref truncate-x i-2))
1707                         len-y low-x-digit))
1708                  (cond ((zerop k) (return))
1709                        (t (decf k)
1710                           (decf low-x-digit)
1711                           (shiftf i i-1 i-2 (1- i-2)))))
1712              q))
1713         ;;; This takes a digit guess, multiplies it by TRUNCATE-Y for a
1714         ;;; result one greater in length than LEN-Y, and subtracts this result
1715         ;;; from TRUNCATE-X. LOW-X-DIGIT is the first digit of X to start
1716         ;;; the subtraction, and we know X is long enough to subtract a LEN-Y
1717         ;;; plus one length bignum from it. Next we check the result of the
1718         ;;; subtraction, and if the high digit in X became negative, then our
1719         ;;; guess was one too big. In this case, return one less than GUESS
1720         ;;; passed in, and add one value of Y back into X to account for
1721         ;;; subtracting one too many. Knuth shows that the guess is wrong on
1722         ;;; the order of 3/b, where b is the base (2 to the digit-size power)
1723         ;;; -- pretty rarely.
1724          (try-bignum-truncate-guess (guess len-y low-x-digit)
1725            (declare (type bignum-index low-x-digit len-y)
1726                     (type bignum-element-type guess))
1727            (let ((carry-digit 0)
1728                  (borrow 1)
1729                  (i low-x-digit))
1730              (declare (type bignum-element-type carry-digit)
1731                       (type bignum-index i)
1732                       (fixnum borrow))
1733              ;; Multiply guess and divisor, subtracting from dividend
1734              ;; simultaneously.
1735              (dotimes (j len-y)
1736                (multiple-value-bind (high-digit low-digit)
1737                    (%multiply-and-add guess
1738                                       (%bignum-ref truncate-y j)
1739                                       carry-digit)
1740                  (declare (type bignum-element-type high-digit low-digit))
1741                  (setf carry-digit high-digit)
1742                  (multiple-value-bind (x temp-borrow)
1743                      (%subtract-with-borrow (%bignum-ref truncate-x i)
1744                                             low-digit
1745                                             borrow)
1746                    (declare (type bignum-element-type x)
1747                             (fixnum temp-borrow))
1748                    (setf (%bignum-ref truncate-x i) x)
1749                    (setf borrow temp-borrow)))
1750                (incf i))
1751              (setf (%bignum-ref truncate-x i)
1752                    (%subtract-with-borrow (%bignum-ref truncate-x i)
1753                                           carry-digit borrow))
1754              ;; See whether guess is off by one, adding one
1755              ;; Y back in if necessary.
1756              (cond ((%digit-0-or-plusp (%bignum-ref truncate-x i))
1757                     guess)
1758                    (t
1759                     ;; If subtraction has negative result, add one
1760                     ;; divisor value back in. The guess was one too
1761                     ;; large in magnitude.
1762                     (let ((i low-x-digit)
1763                           (carry 0))
1764                       (dotimes (j len-y)
1765                         (multiple-value-bind (v k)
1766                             (%add-with-carry (%bignum-ref truncate-y j)
1767                                              (%bignum-ref truncate-x i)
1768                                              carry)
1769                           (declare (type bignum-element-type v))
1770                           (setf (%bignum-ref truncate-x i) v)
1771                           (setf carry k))
1772                         (incf i))
1773                       (setf (%bignum-ref truncate-x i)
1774                             (%add-with-carry (%bignum-ref truncate-x i)
1775                                              0 carry)))
1776                     (%subtract-with-borrow guess 1 1)))))
1777         ;;; This returns the amount to shift y to place a one in the
1778         ;;; second highest bit. Y must be positive. If the last digit
1779         ;;; of y is zero, then y has a one in the previous digit's
1780         ;;; sign bit, so we know it will take one less than digit-size
1781         ;;; to get a one where we want. Otherwise, we count how many
1782         ;;; right shifts it takes to get zero; subtracting this value
1783         ;;; from digit-size tells us how many high zeros there are
1784         ;;; which is one more than the shift amount sought.
1785         ;;;
1786         ;;; Note: This is exactly the same as one less than the
1787         ;;; integer-length of the last digit subtracted from the
1788         ;;; digit-size.
1789         ;;;
1790         ;;; We shift y to make it sufficiently large that doing the
1791         ;;; 2*digit-size by digit-size %FLOOR calls ensures the quotient and
1792         ;;; remainder fit in digit-size.
1793          (shift-y-for-truncate (y)
1794            (let* ((len (%bignum-length y))
1795                   (last (%bignum-ref y (1- len))))
1796              (declare (type bignum-index len)
1797                       (type bignum-element-type last))
1798              (- digit-size (integer-length last) 1)))
1799          ;;; Stores two bignums into the truncation bignum buffers,
1800          ;;; shifting them on the way in. This assumes x and y are
1801          ;;; positive and at least two in length, and it assumes
1802          ;;; truncate-x and truncate-y are one digit longer than x and
1803          ;;; y.
1804          (shift-and-store-truncate-buffers (x len-x y len-y shift)
1805            (declare (type bignum-index len-x len-y)
1806                     (type (integer 0 (#.digit-size)) shift))
1807            (cond ((zerop shift)
1808                   (bignum-replace truncate-x x :end1 len-x)
1809                   (bignum-replace truncate-y y :end1 len-y))
1810                  (t
1811                   (bignum-ashift-left-unaligned x 0 shift (1+ len-x)
1812                                                 truncate-x)
1813                   (bignum-ashift-left-unaligned y 0 shift (1+ len-y)
1814                                                 truncate-y))))) ;; LABELS
1815       ;;; Divide X by Y returning the quotient and remainder. In the
1816       ;;; general case, we shift Y to set up for the algorithm, and we
1817       ;;; use two buffers to save consing intermediate values. X gets
1818       ;;; destructively modified to become the remainder, and we have
1819       ;;; to shift it to account for the initial Y shift. After we
1820       ;;; multiple bind q and r, we first fix up the signs and then
1821       ;;; return the normalized results.
1822       (let* ((x-plusp (%bignum-0-or-plusp x (%bignum-length x)))
1823              (y-plusp (%bignum-0-or-plusp y (%bignum-length y)))
1824              (x (if x-plusp x (negate-bignum x nil)))
1825              (y (if y-plusp y (negate-bignum y nil)))
1826              (len-x (%bignum-length x))
1827              (len-y (%bignum-length y)))
1828         (multiple-value-bind (q r)
1829             (cond ((< len-y 2)
1830                    (bignum-truncate-single-digit x len-x y))
1831                   ((plusp (bignum-compare y x))
1832                    (let ((res (%allocate-bignum len-x)))
1833                      (dotimes (i len-x)
1834                        (setf (%bignum-ref res i) (%bignum-ref x i)))
1835                      (values 0 res)))
1836                   (t
1837                    (let ((len-x+1 (1+ len-x)))
1838                      (setf truncate-x (%allocate-bignum len-x+1))
1839                      (setf truncate-y (%allocate-bignum (1+ len-y)))
1840                      (let ((y-shift (shift-y-for-truncate y)))
1841                        (shift-and-store-truncate-buffers x len-x y
1842                                                          len-y y-shift)
1843                        (values (return-quotient-leaving-remainder len-x+1
1844                                                                   len-y)
1845                                ;; Now that RETURN-QUOTIENT-LEAVING-REMAINDER
1846                                ;; has executed, we just tidy up the remainder
1847                                ;; (in TRUNCATE-X) and return it.
1848                                (cond
1849                                  ((zerop y-shift)
1850                                   (let ((res (%allocate-bignum len-y)))
1851                                     (declare (type bignum-type res))
1852                                     (bignum-replace res truncate-x :end2 len-y)
1853                                     (%normalize-bignum res len-y)))
1854                                  (t
1855                                   (shift-right-unaligned
1856                                    truncate-x 0 y-shift len-y
1857                                    ((= j res-len-1)
1858                                     (setf (%bignum-ref res j)
1859                                           (%ashr (%bignum-ref truncate-x i)
1860                                                  y-shift))
1861                                     (%normalize-bignum res res-len))
1862                                    res))))))))
1863           (let ((quotient (cond ((eq x-plusp y-plusp) q)
1864                                 ((typep q 'fixnum) (the fixnum (- q)))
1865                                 (t (negate-bignum-in-place q))))
1866                 (rem (cond (x-plusp r)
1867                            ((typep r 'fixnum) (the fixnum (- r)))
1868                            (t (negate-bignum-in-place r)))))
1869             (values (if (typep quotient 'fixnum)
1870                         quotient
1871                         (%normalize-bignum quotient (%bignum-length quotient)))
1872                     (if (typep rem 'fixnum)
1873                         rem
1874                         (%normalize-bignum rem (%bignum-length rem))))))))))
1875
1876 \f
1877 ;;;; There used to be a pile of code for implementing division for bignum digits
1878 ;;;; for machines that don't have a 2*digit-size by digit-size divide instruction.
1879 ;;;; This happens to be most machines, but all the SBCL ports seem to be content
1880 ;;;; to implement SB-BIGNUM:%FLOOR as a VOP rather than using the code here.
1881 ;;;; So it's been deleted.  --njf, 2007-02-04
1882 \f
1883 ;;;; general utilities
1884
1885 ;;; Internal in-place operations use this to fixup remaining digits in the
1886 ;;; incoming data, such as in-place shifting. This is basically the same as
1887 ;;; the first form in %NORMALIZE-BIGNUM, but we return the length of the buffer
1888 ;;; instead of shrinking the bignum.
1889 #!-sb-fluid (declaim (sb!ext:maybe-inline %normalize-bignum-buffer))
1890 (defun %normalize-bignum-buffer (result len)
1891   (declare (type bignum-type result)
1892            (type bignum-index len))
1893   (unless (= len 1)
1894     (do ((next-digit (%bignum-ref result (- len 2))
1895                      (%bignum-ref result (- len 2)))
1896          (sign-digit (%bignum-ref result (1- len)) next-digit))
1897         ((not (zerop (logxor sign-digit (%ashr next-digit (1- digit-size))))))
1898         (decf len)
1899         (setf (%bignum-ref result len) 0)
1900         (when (= len 1)
1901               (return))))
1902   len)
1903
1904 ;;; This drops the last digit if it is unnecessary sign information. It repeats
1905 ;;; this as needed, possibly ending with a fixnum. If the resulting length from
1906 ;;; shrinking is one, see whether our one word is a fixnum. Shift the possible
1907 ;;; fixnum bits completely out of the word, and compare this with shifting the
1908 ;;; sign bit all the way through. If the bits are all 1's or 0's in both words,
1909 ;;; then there are just sign bits between the fixnum bits and the sign bit. If
1910 ;;; we do have a fixnum, shift it over for the two low-tag bits.
1911 (defun %normalize-bignum (result len)
1912   (declare (type bignum-type result)
1913            (type bignum-index len)
1914            (inline %normalize-bignum-buffer))
1915   (let ((newlen (%normalize-bignum-buffer result len)))
1916     (declare (type bignum-index newlen))
1917     (unless (= newlen len)
1918       (%bignum-set-length result newlen))
1919     (if (= newlen 1)
1920         (let ((digit (%bignum-ref result 0)))
1921           (if (= (%ashr digit sb!vm:n-positive-fixnum-bits)
1922                  (%ashr digit (1- digit-size)))
1923               (%fixnum-digit-with-correct-sign digit)
1924               result))
1925         result)))
1926
1927 ;;; This drops the last digit if it is unnecessary sign information. It
1928 ;;; repeats this as needed, possibly ending with a fixnum magnitude but never
1929 ;;; returning a fixnum.
1930 (defun %mostly-normalize-bignum (result len)
1931   (declare (type bignum-type result)
1932            (type bignum-index len)
1933            (inline %normalize-bignum-buffer))
1934   (let ((newlen (%normalize-bignum-buffer result len)))
1935     (declare (type bignum-index newlen))
1936     (unless (= newlen len)
1937       (%bignum-set-length result newlen))
1938     result))
1939 \f
1940 ;;;; hashing
1941
1942 ;;; the bignum case of the SXHASH function
1943 (defun sxhash-bignum (x)
1944   (let ((result 316495330))
1945     (declare (type fixnum result))
1946     (dotimes (i (%bignum-length x))
1947       (declare (type index i))
1948       (let ((xi (%bignum-ref x i)))
1949         (mixf result
1950               (logand most-positive-fixnum
1951                       (logxor xi
1952                               (ash xi -7))))))
1953     result))