source: branches/rme-fpe/level-0/l0-bignum64.lisp @ 15779

Last change on this file since 15779 was 13422, checked in by gb, 10 years ago

Simpler BIGNUM-LOGICAL-AND, BIGNUM-LOGICAL-IOR; call directly to the
new LAP functions.
Fix bugs/fenceposts in the ppc64 version of those LAP functions.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 76.4 KB
Line 
1;;;-*- Mode: Lisp; Package: CCL -*-
2;;;
3;;;   Copyright (C) 2009 Clozure Associates
4;;;   Copyright (C) 1994-2001 Digitool, Inc
5;;;   This file is part of Clozure CL. 
6;;;
7;;;   Clozure CL is licensed under the terms of the Lisp Lesser GNU Public
8;;;   License , known as the LLGPL and distributed with Clozure CL as the
9;;;   file "LICENSE".  The LLGPL consists of a preamble and the LGPL,
10;;;   which is distributed with Clozure CL as the file "LGPL".  Where these
11;;;   conflict, the preamble takes precedence. 
12;;;
13;;;   Clozure CL is referenced in the preamble as the "LIBRARY."
14;;;
15;;;   The LLGPL is also available online at
16;;;   http://opensource.franz.com/preamble.html
17
18(in-package "CCL")
19
20#+64-bit-target
21(eval-when (:compile-toplevel :execute)
22  (require "ARCH")
23  (require "NUMBER-MACROS")
24  (require "NUMBER-CASE-MACRO")
25
26  (defsetf bignum-ref bignum-set)
27 
28  (defconstant digit-size 32)
29  (defconstant half-digit-size (/ digit-size 2))
30 
31  (defconstant maximum-bignum-length (1- (ash 1 56)))
32  (defconstant all-ones-digit #xffffffff)
33  (deftype bignum-index () `(integer 0 (,maximum-bignum-length)))
34  (deftype bignum-element-type () `(unsigned-byte ,digit-size))
35  (deftype bignum-half-element-type () `(unsigned-byte ,half-digit-size))
36  (deftype bignum-type () 'bignum)
37  (defmacro %normalize-bignum-macro (big)
38    `(%normalize-bignum-2 t ,big))
39
40  (defmacro %mostly-normalize-bignum-macro (big)
41    `(%normalize-bignum-2 nil ,big))
42  (defmacro %lognot (x)
43    `(logand #xffffffff (lognot (the fixnum ,x))))
44  (defmacro %logior (x y)
45    `(logior (the fixnum ,x) (the fixnum ,y)))
46  (defmacro %logxor (x y)
47    `(logand #xffffffff (logxor (the fixnum ,x) (the fixnum ,y))))
48 
49  ;;; BIGNUM-REPLACE -- Internal.
50  ;;;
51  (defmacro bignum-replace (dest src &key (start1 '0) end1 (start2 '0) end2
52                                 from-end)
53    (once-only ((n-dest dest)
54                (n-src src))
55               (if (and (eq start1 0)(eq start2 0)(null end1)(null end2)(null from-end))
56                 ;; this is all true for some uses today <<
57                 `(%copy-ivector-to-ivector ,n-src 0 ,n-dest 0 (%ilsl 2 (min (the fixnum (%bignum-length ,n-src))
58                                                                         (the fixnum (%bignum-length ,n-dest)))))
59                 (let* ((n-start1 (gensym))
60                        (n-end1 (gensym))
61                        (n-start2 (gensym))
62                        (n-end2 (gensym)))
63                   `(let ((,n-start1 ,start1)
64                          (,n-start2 ,start2)
65                          (,n-end1 ,(or end1 `(%bignum-length ,n-dest)))
66                          (,n-end2 ,(or end2 `(%bignum-length ,n-src))))
67                     ,(if (null from-end)           
68                          `(%copy-ivector-to-ivector
69                            ,n-src (%i* 4 ,n-start2) 
70                            ,n-dest (%i* 4 ,n-start1)
71                            (%i* 4 (min (%i- ,n-end2 ,n-start2) 
72                                    (%i- ,n-end1 ,n-start1))))
73                          `(let ((nwds (min (%i- ,n-end2 ,n-start2)
74                                            (%i- ,n-end1 ,n-start1))))
75                            (%copy-ivector-to-ivector
76                             ,n-src (%ilsl 2 (%i- ,n-end2 nwds))
77                             ,n-dest (%ilsl 2 (%i- ,n-end1 nwds))
78                             (%i* 4 nwds))))))))) 
79 
80
81  ;;;; Shifting.
82 
83  (defconstant all-ones-half-digit #xFFFF) 
84 
85
86;;; %ALLOCATE-BIGNUM must zero all elements.
87;;;
88  (defmacro %allocate-bignum (ndigits)
89    `(%alloc-misc ,ndigits target::subtag-bignum))
90
91  (declaim (inline  %bignum-length))
92
93;;; This macro is used by BIGNUM-ASHIFT-RIGHT,
94;;; BIGNUM-BUFFER-ASHIFT-RIGHT, and BIGNUM-LDB-BIGNUM-RES. They supply
95;;; a termination form that references locals established by this
96;;; form. Source is the source bignum. Start-digit is the first digit
97;;; in source from which we pull bits. Start-pos is the first bit we
98;;; want. Res-len-form is the form that computes the length of the
99;;; resulting bignum. Termination is a DO termination form with a test
100;;; and body. When result is supplied, it is the variable to which
101;;; this binds a newly allocated bignum.
102;;;
103;;; Given start-pos, 1-31 inclusively, of shift, we form the j'th resulting
104;;; digit from high bits of the i'th source digit and the start-pos number of
105;;; bits from the i+1'th source digit.
106  (defmacro shift-right-unaligned (source
107                                   start-digit
108                                   start-pos
109                                   res-len-form
110                                   termination
111                                   &optional result)
112    `(let* ((high-bits-in-first-digit (- digit-size ,start-pos))
113            (res-len ,res-len-form)
114            (res-len-1 (1- res-len))
115            ,@(if result `((,result (%allocate-bignum res-len)))))
116      (declare (type bignum-index res-len res-len-1))
117      (do ((i ,start-digit i+1)
118           (i+1 (1+ ,start-digit) (1+ i+1))
119           (j 0 (1+ j)))
120          ,termination
121        (declare (type bignum-index i i+1 j))
122        (setf (bignum-ref ,(if result result source) j)
123              (%logior (%digit-logical-shift-right (bignum-ref ,source i)
124                                                   ,start-pos)
125                       (%ashl (bignum-ref ,source i+1)
126                              high-bits-in-first-digit))))))
127
128
129  )
130
131
132#+64-bit-target
133(progn
134
135;;; Extract the length of the bignum.
136;;;
137(defun %bignum-length (bignum)
138  (uvsize bignum)) 
139
140
141
142;;; We can probably do better than UVREF here, but
143;;; a) it's not -that- bad
144;;; b) it does some bounds/sanity checking, which isn't a bad idea.
145
146(defmacro bignum-ref (b i)
147  `(%typed-miscref :bignum ,b ,i))
148
149(defmacro bignum-set (b i val)
150  `(%typed-miscset :bignum ,b ,i ,val))
151
152
153(defun bignum-plusp (b)
154  (not (logbitp (1- digit-size) (the bignum-element-type (bignum-ref b (1- (%bignum-length b)))))))
155
156;;; Return T if digit is positive, or NIL if negative.
157(defun %digit-0-or-plusp (digit)
158  (declare (type bignum-element-type digit))
159  (not (logbitp (1- digit-size) digit)))
160
161(defun %bignum-0-or-plusp (bignum len)
162  (declare (type bignum-type bignum)
163           (type bignum-index len))
164  (%digit-0-or-plusp (bignum-ref bignum (1- len))))
165
166(defun bignum-minusp (b)
167  (logbitp 31 (the fixnum (bignum-ref b (1- (%bignum-length b))))))
168
169(defun %sign-digit (b i)
170  (%ashr (bignum-ref b (1- i)) (1- digit-size)))
171
172;;; Return the sign of bignum (0 or -1) as a fixnum
173(defun %bignum-sign (b)
174  (if (logbitp 31 (the fixnum (bignum-ref b (1- (%bignum-length b)))))
175    -1
176    0))
177
178         
179(defun %add-with-carry (a-digit b-digit carry-in)
180  (declare (fixnum a-digit b-digit carry-in))
181  (setq a-digit (logand all-ones-digit a-digit)
182        b-digit (logand all-ones-digit b-digit))
183  (let* ((sum (+ carry-in (the fixnum (+ a-digit b-digit)))))
184    (declare (fixnum sum))
185    (values (logand all-ones-digit sum) (logand 1 (the fixnum (ash sum -32))))))
186
187(defun %subtract-with-borrow (a-digit b-digit borrow-in)
188  (declare (fixnum a-digit b-digit borrow-in))
189  (setq a-digit (logand all-ones-digit a-digit)
190        b-digit (logand all-ones-digit b-digit))
191  (let* ((diff (- (the fixnum (- a-digit b-digit))
192                  (the fixnum (- 1 borrow-in)))))
193    (declare (fixnum diff))
194    (values (logand all-ones-digit diff)
195            (- 1 (logand (the fixnum (ash diff -32)) 1)))))
196
197
198
199(defun %compare-digits (bignum-a bignum-b idx)
200  (let* ((a (bignum-ref bignum-a idx))
201         (b (bignum-ref bignum-b idx)))
202    (declare (fixnum a b))
203    (if (= a b)
204      0
205      (if (> a b)
206        1
207        -1))))
208
209
210;;;; Addition.
211(defun add-bignums (a b)
212  (let* ((len-a (%bignum-length a))
213         (len-b (%bignum-length b)))
214    (declare (bignum-index len-a len-b)
215             (optimize (speed 3) (safety 0)))
216    (when (> len-b len-a)
217      (rotatef a b)
218      (rotatef len-a len-b))
219    (let* ((len-res (1+ len-a))
220           (res (%allocate-bignum len-res))
221           (carry 0)
222           (sign-b (%bignum-sign b)))
223        (dotimes (i len-b)
224          (let* ((sum (+
225                       (the fixnum (+ (the bignum-element-type (bignum-ref a i))
226                                      (the bignum-element-type (bignum-ref b i))))
227                       carry)))
228            (declare (fixnum sum))
229            (setf (bignum-ref res i) sum)
230            (setq carry (logand 1 (the fixnum (ash sum -32))))))
231        (if (/= len-a len-b)
232          (finish-bignum-add  res carry a sign-b len-b len-a)
233          (setf (bignum-ref res len-a)
234                (+ (the fixnum carry)
235                   (the fixnum (+ (the bignum-element-type (%bignum-sign a))
236                                  sign-b)))))
237        (%normalize-bignum-macro res))))
238
239(defun add-bignum-and-fixnum (bignum fixnum)
240  (declare (bignum-type bignum)
241           (fixnum fixnum)
242           (optimize (speed 3) (safety 0)))
243  (let* ((len-bignum (%bignum-length bignum))
244         (len-res (1+ len-bignum))
245         (res (%allocate-bignum len-res))
246         (low (logand all-ones-digit fixnum))
247         (high (logand all-ones-digit (the fixnum (ash fixnum -32)))))
248    (declare (bignum-index len-bignum)
249             (bignum-type res)
250             (bignum-element-type low high))
251    (let* ((sum0 (+ (the bignum-element-type (bignum-ref bignum 0)) low))
252           (sum1 (+ (the fixnum (+ (the bignum-element-type (bignum-ref bignum 1))
253                                   high))
254                    (the fixnum (logand 1 (ash sum0 -32)))))
255           (carry (logand 1 (ash sum1 -32))))
256      (declare (fixnum sum0 sum1) (bignum-element-type carry))
257      (setf (bignum-ref res 0) sum0
258            (bignum-ref res 1) sum1)
259      (if (> len-bignum 2)
260        (finish-bignum-add  res carry bignum (ash fixnum (- (- target::nbits-in-word target::fixnumshift))) 2 len-bignum)
261        (setf (bignum-ref res 2)
262              (+ (the fixnum carry)
263                 (the fixnum (+ (the bignum-element-type (%bignum-sign bignum))
264                                (the fixnum (ash fixnum (- (- target::nbits-in-word target::fixnumshift)))))))))
265      (%normalize-bignum-macro res))))
266
267
268
269
270
271;;; B was shorter than A; keep adding B's sign digit to each remaining
272;;; digit of A, propagating the carry.
273(defun finish-bignum-add (result carry a sign-b start end)
274  (declare (type bignum-index start end)
275           (bignum-element-type sign-b carry)
276           (optimize (speed 3) (safety 0)))
277  (do* ((i start (1+ i))
278        (sign-b (logand all-ones-digit sign-b)))
279       ((= i end)
280        (setf (bignum-ref result end)
281              (the fixnum (+
282                           (the fixnum (+ (the fixnum
283                                            (logand all-ones-digit
284                                                    (the fixnum
285                                                      (%sign-digit a end))))
286                                          sign-b))
287                           carry))))
288    (declare (fixnum i) (bignum-element-type sign-b))
289    (let* ((sum (the fixnum (+ (the fixnum (+ (bignum-ref a i)
290                                              sign-b))
291                               carry))))
292      (declare (fixnum sum))
293      (setf (bignum-ref result i) sum)
294      (setq carry (logand 1 (the fixnum (ash sum -32)))))))
295
296
297
298
299;;;; Subtraction.
300(defun subtract-bignum (a b)
301  (let* ((len-a (%bignum-length a))
302         (len-b (%bignum-length b))
303         (len-res (1+ (max len-a len-b)))
304         (res (%allocate-bignum len-res)))
305    (declare (bignum-index len-a len-b len-res))
306    (bignum-subtract-loop a len-a b len-b res)
307    (%normalize-bignum-macro res)))
308
309(defun bignum-subtract-loop (a len-a b len-b res)
310  (declare (bignum-index len-a len-b )
311           (optimize (speed 3) (safety 0)))
312  (let* ((len-res (%bignum-length res)))
313    (declare (bignum-index len-res))
314    (let* ((borrow 1)
315           (sign-a (%bignum-sign a))
316           (sign-b (%bignum-sign b)))
317      (declare (bignum-element-type borrow sign-a sign-b))
318      (dotimes (i (the bignum-index len-res))
319        (multiple-value-bind (result-digit borrow-out)
320            (%subtract-with-borrow
321             (if (< i len-a)
322               (bignum-ref a i)
323               sign-a)
324             (if (< i len-b)
325               (bignum-ref b i)
326               sign-b)
327             borrow)
328          (setf (bignum-ref res i) result-digit
329                borrow borrow-out))))))
330
331
332;;;; Multiplication.
333
334#||
335;;; These parameters match GMP's.
336(defvar *sqr-basecase-threshold* 5)
337(defvar *sqr-karatsuba-threshold* 22)
338(defvar *mul-karatsuba-threshold* 10)
339
340;;; Squaring is often simpler than multiplication.  This should never
341;;; be called with (>= N *sqr-karatsuba-threshold*).
342(defun mpn-sqr-basecase (prodp up n)
343  (declare (fixnum prodp up n))
344  (declare (optimize (speed 3) (safety 0) (space 0)))
345  (umulppm up up prodp)
346  (when (> n 1)
347    (%stack-block ((tarr (* 4 (* 2 *sqr-karatsuba-threshold*))))
348      (let* ((tp (macptr->fixnum tarr)))
349        (mpn-mul-1 tp
350                   (the fixnum (1+ up))
351                   (the fixnum (1- n))
352                   up
353                   (the fixnum (+ tp (the fixnum (1- n)))))
354        (do* ((i 2 (1+ i)))
355             ((= i n))
356          (declare (fixnum i))
357          (mpn-addmul-1 (the fixnum (- (the fixnum (+ tp (the fixnum (+ i i))))
358                                       2))
359                        (the fixnum (+ up i))
360                        (the fixnum (- n i))
361                        (the fixnum (+ up (the fixnum (1- i))))
362                        (the fixnum (+ tp (the fixnum (+ n (the fixnum (- i 2))))))))
363        (do* ((i 1 (1+ i))
364              (ul (1+ up) (1+ ul)))
365             ((= i n))
366          (declare (fixnum i ul))
367          (umulppm ul ul (the fixnum (+ prodp (the fixnum (+ i i))))))
368        (let* ((2n-2 (- (the fixnum (+ n n)) 2))
369               (carry (mpn-lshift-1 tp tp 2n-2)))
370          (declare (fixnum 2n-2 carry))
371          (incf carry (the fixnum (mpn-add-n (the fixnum (1+ prodp))
372                                             (the fixnum (1+ prodp))
373                                             tp
374                                             2n-2)))
375          (add-fixnum-to-limb carry (the fixnum (+ prodp
376                                                   (the fixnum (1-
377                                                                (the fixnum
378                                                                  (+ n n))))))))))))
379
380;;; For large enough values of N, squaring via Karatsuba-style
381;;; divide&conquer is faster than in the base case.
382(defun mpn-kara-sqr-n (p a n ws)
383  (declare (fixnum p a n ws))
384  (declare (optimize (speed 3) (safety 0) (space 0)))
385  (%stack-block ((limbs 16))
386    (let* ((w (macptr->fixnum limbs))
387           (w0 (1+ w))
388           (w1 (1+ w0))
389           (xx (1+ w1))
390           (n2 (ash n -1))
391           (x 0)
392           (y 0)
393           (i 0))
394      (declare (fixnum w w0 w1 xx n2 x y i))
395      (cond ((logbitp 0 n)
396             ;; Odd length
397             (let* ((n3 (- n n2))
398                    (n1 0)
399                    (nm1 0))
400               (declare (fixnum n3 n1 nm1))
401               (copy-limb (the fixnum (+ a n2)) w)
402               (if (not (limb-zerop w))
403                 (add-fixnum-to-limb
404                  (the fixnum
405                    (- (the fixnum (mpn-sub-n p a (the fixnum (+ a n3)) n2))))
406                  w)
407                 (progn
408                   (setq i n2)
409                   (loop
410                     (decf i)
411                     (copy-limb (the fixnum (+ a i)) w0)
412                     (copy-limb (the fixnum (+ a (the fixnum (+ n3 i)))) w1)
413                     (if (or (not (zerop (the fixnum (compare-limbs w0 w1))))
414                             (= i 0))
415                       (return)))
416                   (if (< (the fixnum (compare-limbs w0 w1)) 0)
417                     (setq x (+ a n3)
418                           y a)
419                     (setq y (+ a n3)
420                           x a))
421                   (mpn-sub-n p x y n2)))
422               (copy-limb w (the fixnum (+ p n2)))
423               (setq n1 (1+ n))
424               (cond ((< n3 *sqr-basecase-threshold*)
425                      (mpn-mul-basecase ws p n3 p n3)
426                      (mpn-mul-basecase p a n3 a n3))
427                     ((< n3 *sqr-karatsuba-threshold*)
428                      (mpn-sqr-basecase ws p n3)
429                      (mpn-sqr-basecase p a n3))
430                     (t
431                      (mpn-kara-sqr-n ws p n3 (the fixnum (+ ws n1)))
432                      (mpn-kara-sqr-n p  a n3 (the fixnum (+ ws n1)))))
433               (cond ((< n2 *sqr-basecase-threshold*)
434                      (mpn-mul-basecase (the fixnum (+ p n1))
435                                        (the fixnum (+ a n3))
436                                        n2
437                                        (the fixnum (+ a n3))
438                                        n2))
439                     ((< n2 *sqr-karatsuba-threshold*)
440                      (mpn-sqr-basecase (the fixnum (+ p n1))
441                                        (the fixnum (+ a n3))
442                                        n2))
443                     (t
444                      (mpn-kara-sqr-n (the fixnum (+ p n1))
445                                      (the fixnum (+ a n3))
446                                      n2
447                                      (the fixnum (+ ws n1)))))
448               (mpn-sub-n ws p ws n1)
449               (setq nm1 (1- n))
450               (unless (zerop (the fixnum
451                                (mpn-add-n ws
452                                           (the fixnum (+ p n1))
453                                           ws
454                                           nm1)))
455                 (copy-limb (the fixnum (+ ws nm1)) xx)
456                 (add-fixnum-to-limb 1 xx)
457                 (copy-limb xx (the fixnum (+ ws nm1)))
458                 (if (limb-zerop xx)
459                   (add-fixnum-to-limb 1 (the fixnum (+ ws n)))))
460               (unless (zerop
461                        (the fixnum
462                          (mpn-add-n (the fixnum (+ p n3))
463                                     (the fixnum (+ p n3))
464                                     ws
465                                     n1)))
466                 (mpn-incr-u (the fixnum (+ p (the fixnum (+ n1 n3))))
467                             1))))
468            (t ; N is even
469             (setq i n2)
470             (loop
471               (decf i)
472               (copy-limb (the fixnum (+ a i)) w0)
473               (copy-limb (the fixnum (+ a (the fixnum (+ n2 i)))) w1)
474               (if (or (not (zerop (the fixnum (compare-limbs w0 w1))))
475                       (= i 0))
476                 (return)))
477             (if (< (the fixnum (compare-limbs w0 w1)) 0)
478               (setq x (+ a n2)
479                     y a)
480               (setq y (+ a n2)
481                     x a))
482             (mpn-sub-n p x y n2)
483             (cond ((< n2 *sqr-basecase-threshold*)
484                    (mpn-mul-basecase ws p n2 p n2)
485                    (mpn-mul-basecase p a n2 a n2)
486                    (mpn-mul-basecase (the fixnum (+ p n))
487                                      (the fixnum (+ a n2))
488                                      n2
489                                      (the fixnum (+ a n2))
490                                      n2))
491                   ((< n2 *sqr-karatsuba-threshold*)
492                    (mpn-sqr-basecase ws p n2)
493                    (mpn-sqr-basecase p a n2)
494                    (mpn-sqr-basecase (the fixnum (+ p n))
495                                      (the fixnum (+ a n2))
496                                      n2))
497                   (t
498                    (mpn-kara-sqr-n ws p n2 (the fixnum (+ ws n)))
499                    (mpn-kara-sqr-n p  a n2 (the fixnum (+ ws n)))
500                    (mpn-kara-sqr-n (the fixnum (+ p n))
501                                    (the fixnum (+ a n2))
502                                    n2
503                                    (the fixnum (+ ws n)))))
504             (let* ((ww (- (the fixnum (mpn-sub-n ws p ws n)))))
505               (declare (fixnum ww))
506               (setq ww (+ ww (mpn-add-n ws (the fixnum (+ p n)) ws n)))
507               (setq ww (+ ww (mpn-add-n (the fixnum (+ p n2))
508                                         (the fixnum (+ p n2))
509                                         ws
510                                         n)))
511               (mpn-incr-u (the fixnum (+ p (the fixnum (+ n2 n)))) ww)))))))
512
513;;; Karatsuba subroutine: multiply A and B, store result at P, use WS
514;;; as scrach space.  Treats A and B as if they were both of size N;
515;;; if that's not true, caller must fuss around the edges.
516(defun mpn-kara-mul-n (p a b n ws)
517  (declare (fixnum p a b n ws))
518  (declare (optimize (speed 3) (safety 0) (space 0)))
519  (%stack-block ((limbs 16))
520    (let* ((w (macptr->fixnum limbs))
521           (w0 (1+ w))
522           (w1 (1+ w0))
523           (xx (1+ w1))
524           (x 0)
525           (y 0)
526           (i 0)
527           (n2 (ash n -1))
528           (sign 0))
529      (declare (fixnum w w0 w1 xx x y i n2 sign))
530      (cond ((logbitp 0 n)
531             (let* ((n1 0)
532                    (n3 (- n n2))
533                    (nm1 0))
534               (declare (fixnum n1 n3 nm1))
535               (copy-limb (the fixnum (+ a n2)) w)
536               (if (not (limb-zerop w))
537                 (add-fixnum-to-limb
538                  (the fixnum (- (mpn-sub-n p a (the fixnum (+ a n3)) n2))) w)
539                 (progn
540                   (setq i n2)
541                   (loop
542                     (decf i)
543                     (copy-limb (the fixnum (+ a i)) w0)
544                     (copy-limb (the fixnum (+ a (the fixnum (+ n3 i)))) w1)
545                     (if (or (not (zerop (the fixnum (compare-limbs w0 w1))))
546                             (zerop i))
547                       (return)))
548                   (if (< (the fixnum (compare-limbs w0 w1)) 0)
549                     (setq x (+ a n3)
550                           y a
551                           sign -1)
552                     (setq x a
553                           y (+ a n3)))
554                   (mpn-sub-n p x y n2)))
555               (copy-limb w (the fixnum (+ p n2)))
556               (copy-limb (the fixnum (+ b n2)) w)
557               (if (not (limb-zerop w))
558                 (add-fixnum-to-limb
559                  (the fixnum (- (the fixnum (mpn-sub-n (the fixnum (+ p n3))
560                                                        b
561                                                        (the fixnum (+ b n3))
562                                                        n2))))
563                  w)
564                 (progn
565                   (setq i n2)
566                   (loop
567                     (decf i)
568                     (copy-limb (the fixnum (+ b i)) w0)
569                     (copy-limb (the fixnum (+ b (the fixnum (+ n3 i)))) w1)
570                     (if (or (not (zerop (the fixnum (compare-limbs w0 w1))))
571                             (zerop i))
572                       (return)))
573                   (if (< (the fixnum (compare-limbs w0 w1)) 0)
574                     (setq x (+ b n3)
575                           y b
576                           sign (lognot sign))
577                     (setq x b
578                           y (+ b n3)))
579                   (mpn-sub-n (the fixnum (+ p n3)) x y n2)))
580               (copy-limb w (the fixnum (+ p n)))
581               (setq n1 (1+ n))
582               (cond
583                 ((< n2 *mul-karatsuba-threshold*)
584                  (cond
585                    ((< n3 *mul-karatsuba-threshold*)
586                     (mpn-mul-basecase ws p n3 (the fixnum (+ p n3)) n3)
587                     (mpn-mul-basecase p a n3 b n3))
588                    (t
589                     (mpn-kara-mul-n ws p (the fixnum (+ p n3)) n3 (the fixnum (+ ws n1)))
590                     (mpn-kara-mul-n p a b n3 (the fixnum (+ ws n1)))))
591                  (mpn-mul-basecase (the fixnum (+ p n1))
592                                    (the fixnum (+ a n3))
593                                    n2
594                                    (the fixnum (+ b n3))
595                                    n2))
596                 (t
597                  (mpn-kara-mul-n ws p (the fixnum (+ p n3)) n3 (the fixnum (+ ws n1)))
598                  (mpn-kara-mul-n p a b n3 (the fixnum (+ ws n1)))
599                  (mpn-kara-mul-n (the fixnum (+ p n1))
600                                  (the fixnum (+ a n3))
601                                  (the fixnum (+ b n3))
602                                  n2
603                                  (the fixnum (+ ws n1)))))
604               (if (not (zerop sign))
605                 (mpn-add-n ws p ws n1)
606                 (mpn-sub-n ws p ws n1))
607               (setq nm1 (1- n))
608               (unless (zerop (the fixnum (mpn-add-n ws
609                                                     (the fixnum (+ p n1))
610                                                     ws
611                                                     nm1)))
612                 (copy-limb (the fixnum (+ ws nm1)) xx)
613                 (add-fixnum-to-limb 1 xx)
614                 (copy-limb xx (the fixnum (+ ws nm1)))
615                 (if (limb-zerop xx)
616                   (add-fixnum-to-limb 1 (the fixnum (+ ws n)))))
617               (unless (zerop (the fixnum
618                                (mpn-add-n (the fixnum (+ p n3))
619                                           (the fixnum (+ p n3))
620                                           ws
621                                           n1)))
622                 (mpn-incr-u (the fixnum
623                               (+ p (the fixnum (+ n1 n3)))) 1))))
624            (t                          ; even length
625             (setq i n2)
626             (loop
627               (decf i)
628               (copy-limb (the fixnum (+ a i)) w0)
629               (copy-limb (the fixnum (+ a (the fixnum (+ n2 i)))) w1)
630               (if (or (not (zerop (the fixnum (compare-limbs w0 w1))))
631                       (zerop i))
632                 (return)))
633             (setq sign 0)
634             (if (< (the fixnum (compare-limbs w0 w1)) 0)
635               (setq x (+ a n2)
636                     y a
637                     sign -1)
638               (setq x a
639                     y (+ a n2)))
640             (mpn-sub-n p x y n2)
641             (setq i n2)
642             (loop
643               (decf i)
644               (copy-limb (the fixnum (+ b i)) w0)
645               (copy-limb (the fixnum (+ b (the fixnum (+ n2 i)))) w1)
646               (if (or (not (zerop (the fixnum (compare-limbs w0 w1))))
647                       (zerop i))
648                 (return)))           
649             (if (< (the fixnum (compare-limbs w0 w1)) 0)
650               (setq x (+ b n2)
651                     y b
652                     sign (lognot sign))
653               (setq x b
654                     y (+ b n2)))
655             (mpn-sub-n (the fixnum (+ p n2)) x y n2)
656             (cond
657               ((< n2 *mul-karatsuba-threshold*)
658                (mpn-mul-basecase ws p n2 (the fixnum (+ p n2)) n2)
659                (mpn-mul-basecase p a n2 b n2)
660                (mpn-mul-basecase (the fixnum (+ p n))
661                                  (the fixnum (+ a n2))
662                                  n2
663                                  (the fixnum (+ b n2))
664                                  n2))
665               (t
666                (mpn-kara-mul-n ws p (the fixnum (+ p n2)) n2
667                                (the fixnum (+ ws n)))
668                (mpn-kara-mul-n p a b n2 (the fixnum (+ ws n)))
669                (mpn-kara-mul-n (the fixnum (+ p n))
670                                (the fixnum (+ a n2))
671                                (the fixnum (+ b n2))
672                                n2
673                                (the fixnum (+ ws n)))))
674             (let* ((ww (if (not (zerop sign))
675                          (mpn-add-n ws p ws n)
676                          (- (the fixnum (mpn-sub-n ws p ws n))))))
677               (declare (fixnum ww))
678               (setq ww (+ ww (mpn-add-n ws (the fixnum (+ p n)) ws n)))
679               (setq ww (+ ww (mpn-add-n (the fixnum (+ p n2))
680                                         (the fixnum (+ p n2))
681                                         ws
682                                         n)))
683               (mpn-incr-u (the fixnum (+ p (the fixnum (+ n2 n)))) ww)))))))
684
685;;; Square UP, of length UN.  I wonder if a Karatsuba multiply might be
686;;; faster than a basecase square.
687(defun mpn-sqr-n (prodp up un)
688  (declare (fixnum prodp up un))
689  (declare (optimize (speed 3) (safety 0) (space 0)))
690  (if (< un *sqr-basecase-threshold*)
691    (mpn-mul-basecase prodp up un up un)
692    (if (< un *sqr-karatsuba-threshold*)
693      (mpn-sqr-basecase prodp up un)
694      (%stack-block ((wsptr (mpn-kara-sqr-n-tsize un)))
695        (mpn-kara-sqr-n prodp up un (macptr->fixnum wsptr))))))
696
697;;; Subroutine: store AxB at P.  Assumes A & B to be of length N
698(defun mpn-mul-n (p a b n)
699  (declare (fixnum p a b n))
700  (declare (optimize (speed 3) (safety 0) (space 0))) 
701  (if (< n *mul-karatsuba-threshold*)
702    (mpn-mul-basecase p a n b n)
703    (%stack-block ((wsptr (mpn-kara-mul-n-tsize n)))
704      (mpn-kara-mul-n p a b n (macptr->fixnum wsptr)))))
705
706
707;;; Multiply [UP,UN] by [VP,VN].  UN must not be less than VN.
708;;; This does Karatsuba if operands are big enough; if they are
709;;; and they differ in size, this computes the product of the
710;;; smaller-size slices, then fixes up the resut.
711(defun mpn-mul (prodp up un vp vn)
712  (declare (fixnum prodp up un vp vn))
713  (declare (optimize (speed 3) (safety 0) (space 0)))
714  ;(assert (>= un vn 1))
715  (if (and (= up vp) (= un vn))
716    (mpn-sqr-n prodp up un)
717    (if (< vn *mul-karatsuba-threshold*)
718      (mpn-mul-basecase prodp up un vp vn)
719      (let* ((l vn))
720        (declare (fixnum l))
721        (mpn-mul-n prodp up vp vn)
722        (unless (= un vn)
723          (incf prodp vn)
724          (incf up vn)
725          (decf un vn)
726          (if (< un vn)
727            (psetq un vn vn un up vp vp up))
728          (%stack-block ((wsptr
729                          (the fixnum
730                            (+ 8
731                               (the fixnum
732                                 (* 4
733                                    (the fixnum
734                                      (+ vn
735                                         (if (>= vn *mul-karatsuba-threshold*)
736                                           vn
737                                           un)))))))))
738            (setf (%get-unsigned-long wsptr 0) 0
739                  (%get-unsigned-long wsptr 4) 0)
740            (let* ((tt (macptr->fixnum wsptr))
741                   (c (1+ tt))
742                   (ws (1+ c)))
743              (declare (fixnum tt c ws ))
744              (do* ()
745                   ((< vn *mul-karatsuba-threshold*))
746                (mpn-mul-n ws up vp vn)
747                (cond ((<= l (the fixnum (+ vn vn)))
748                       (add-fixnum-to-limb (mpn-add-n prodp prodp ws l) tt)
749                       (unless (= l (the fixnum (+ vn vn)))
750                         (copy-fixnum-to-limb
751                          (mpn-add-1 (the fixnum (+ prodp l))
752                                     (the fixnum (+ ws l))
753                                     (the fixnum (- (the fixnum (+ vn vn)) l))
754                                     tt)
755                          tt)
756                         (setq l (the fixnum (+ vn vn)))))
757                      (t
758                       (copy-fixnum-to-limb
759                        (mpn-add-n prodp prodp ws (the fixnum (+ vn vn))) c)
760                       (add-fixnum-to-limb
761                        (mpn-add-1 (the fixnum (+ prodp (the fixnum (+ vn vn))))
762                                   (the fixnum (+ prodp (the fixnum (+ vn vn))))
763                                   (the fixnum (- l (the fixnum (+ vn vn))))
764                                   c)
765                        tt)))
766                (incf prodp vn)
767                (decf l vn)
768                (incf up vn)
769                (decf un vn)
770                (if (< un vn)
771                  (psetq up vp vp up un vn vn un)))
772              (unless (zerop vn)
773                (mpn-mul-basecase ws up un vp vn)
774                (cond ((<= l (the fixnum (+ un vn)))
775                       (add-fixnum-to-limb
776                        (mpn-add-n prodp prodp ws l)
777                        tt)
778                       (unless (= l (the fixnum (+ un vn)))
779                         (copy-fixnum-to-limb
780                          (mpn-add-1 (the fixnum (+ prodp l))
781                                     (the fixnum (+ ws l))
782                                     (the fixnum (- (the fixnum (+ un vn)) l))
783                                     tt)
784                          tt)))
785                      (t
786                       (copy-fixnum-to-limb
787                        (mpn-add-n prodp prodp ws (the fixnum (+ un vn)))
788                        c)
789                       (add-fixnum-to-limb
790                        (mpn-add-1
791                         (the fixnum (+ prodp (the fixnum (+ un vn))))
792                         (the fixnum (+ prodp (the fixnum (+ un vn))))
793                         (the fixnum (- (the fixnum (- l un)) vn))
794                         c)
795                        tt)))))))))))
796||#
797
798(defun multiply-bignums (a b)
799  (let* ((signs-differ (not (eq (bignum-minusp a) (bignum-minusp b)))))
800    (flet ((multiply-unsigned-bignums (a b)
801             (let* ((len-a (%bignum-length a))
802                    (len-b (%bignum-length b))
803                    (len-res (+ len-a len-b))
804                    (res (%allocate-bignum len-res)))
805               (declare (bignum-index len-a len-b len-res))
806               (dotimes (i len-a)
807                 (declare (type bignum-index i))
808                 (%multiply-and-add-loop a b res i len-b))
809               res))
810           (multiply-unsigned-bignums64 (a b)
811             (let* ((len-a (ceiling (%bignum-length a) 2))
812                    (len-b (ceiling (%bignum-length b) 2))
813                    (len-res (+ len-a len-b))
814                    (res (%allocate-bignum (+ len-res len-res))))
815               (declare (bignum-index len-a len-b len-res))
816               (dotimes (i len-a)
817                 (declare (type bignum-index i))
818                 (%multiply-and-add-loop64 a b res i len-b))
819               res)))
820      (let* ((res (with-negated-bignum-buffers a b
821                                               multiply-unsigned-bignums64)))
822        (if signs-differ (negate-bignum-in-place res))
823        (%normalize-bignum-macro res)))))
824
825(defun multiply-bignum-and-fixnum (bignum fixnum)
826  (declare (type bignum-type bignum) (fixnum fixnum))
827  (if (eql fixnum 1)
828    bignum
829    (with-small-bignum-buffers ((big-fix fixnum))
830      (multiply-bignums bignum big-fix))))
831
832#+slower
833(defun multiply-bignum-and-fixnum (bignum fixnum)
834  (declare (type bignum-type bignum) (fixnum fixnum))
835  (if (eql fixnum 1)
836    bignum
837    (if (eql fixnum target::target-most-negative-fixnum)
838      (with-small-bignum-buffers ((big-fix fixnum))
839        (multiply-bignums bignum big-fix))
840      (let* ((big-len (%bignum-length bignum))
841             (big-neg (bignum-minusp bignum))
842             (signs-differ (not (eq big-neg (minusp fixnum)))))
843        (flet ((multiply-unsigned-bignum-and-2-digit-fixnum (a len-a high low)
844                 (declare (bignum-type a)
845                          (bignum-element-type high low)
846                          (bignum-index len-a)
847                          (optimize (speed 3) (safety 0)))
848                 (let* ((len-res (+ len-a 2))
849                        (res (%allocate-bignum len-res)) )
850                   (declare (bignum-index len-a  len-res))
851                   (dotimes (i len-a)
852                     (declare (type bignum-index i))
853                     (let* ((carry-digit 0)
854                            (x (bignum-ref a i))
855                            (k i))
856                       (declare (fixnum k))
857                       (multiple-value-bind (big-carry res-digit)
858                           (%multiply-and-add4 x
859                                               low
860                                               (bignum-ref res k)
861                                               carry-digit)
862                         (setf (bignum-ref res k) res-digit
863                               carry-digit big-carry
864                               k (1+ k)))
865                       (multiple-value-bind (big-carry res-digit)
866                           (%multiply-and-add4 x
867                                               high
868                                               (bignum-ref res k)
869                                               carry-digit)
870                         (setf (bignum-ref res k) res-digit
871                               carry-digit big-carry
872                               k (1+ k)))
873                       (setf (bignum-ref res k) carry-digit)))
874                   res))
875               (multiply-unsigned-bignum-and-1-digit-fixnum (a len-a fix)
876                 (declare (bignum-type a)
877                          (bignum-element-type fix)
878                          (bignum-index len-a)
879                          (optimize (speed 3) (safety 0)))
880                 (let* ((len-res (+ len-a 1))
881                        (res (%allocate-bignum len-res)) )
882                   (declare (bignum-index len-a  len-res))
883                   (dotimes (i len-a)
884                     (declare (type bignum-index i))
885                     (let* ((carry-digit 0)
886                            (x (bignum-ref a i))
887                            (k i))
888                       (declare (fixnum k))
889                       (multiple-value-bind (big-carry res-digit)
890                           (%multiply-and-add4 x
891                                               fix
892                                               (bignum-ref res k)
893                                               carry-digit)
894                         (setf (bignum-ref res k) res-digit
895                               carry-digit big-carry
896                               k (1+ k)))
897                       (setf (bignum-ref res k) carry-digit)))
898                   res)))
899          (let* ((fixnum (if (< fixnum 0) (- fixnum) fixnum))
900                 (low (logand (1- (ash 1 32)) fixnum))
901                 (high (unless (<= (%fixnum-intlen fixnum) 32)
902                         (ldb (byte 32 32) fixnum)))
903                 (res (if big-neg
904                        (let* ((neg-len (1+ big-len)))
905                          (declare (type bignum-index neg-len))
906                          (with-bignum-buffers ((neg neg-len))
907                            (negate-bignum bignum nil neg)
908                            (setq neg-len (%bignum-length bignum))
909                            (if high
910                              (multiply-unsigned-bignum-and-2-digit-fixnum
911                               neg
912                               neg-len
913                               high
914                               low)
915                              (multiply-unsigned-bignum-and-1-digit-fixnum
916                               neg
917                               neg-len
918                               low))))
919                        (if high
920                          (multiply-unsigned-bignum-and-2-digit-fixnum
921                           bignum
922                           big-len
923                           high
924                           low)
925                          (multiply-unsigned-bignum-and-1-digit-fixnum
926                           bignum
927                           big-len
928                           low)))))
929            (if signs-differ (negate-bignum-in-place res))
930            (%normalize-bignum-macro res)))))))
931
932
933;; assume we already know result won't fit in a fixnum
934;; only caller is fixnum-*-2
935;;
936
937(defun multiply-fixnums (a b)
938  (declare (fixnum a b))
939  (* a b))
940
941
942;;;; GCD.
943
944
945;;; Both args are > 0.
946(defun bignum-fixnum-gcd (bignum fixnum)
947  (let* ((rem (bignum-truncate-by-fixnum-no-quo bignum fixnum)))
948    (declare (fixnum rem))
949    (if (zerop rem)
950      fixnum
951      (%fixnum-gcd rem fixnum))))
952
953
954
955;;; NEGATE-BIGNUM -- Public.
956;;;
957;;; Fully-normalize is an internal optional.  It cause this to always return
958;;; a bignum, without any extraneous digits, and it never returns a fixnum.
959;;;
960(defun negate-bignum (x &optional (fully-normalize t) res)
961  (declare (type bignum-type x))
962  (let* ((len-x (%bignum-length x))
963         (len-res (1+ len-x))
964         (minusp (bignum-minusp x))
965         (res (or res (%allocate-bignum len-res))))
966    (declare (type bignum-index len-x len-res)) ;Test len-res for range?
967    (let ((carry (bignum-negate-loop-really x len-x res)))
968      (declare (fixnum carry))
969      (if (zerop carry)
970        (setf (bignum-ref res len-x) (if minusp 0 all-ones-digit))
971        (setf (bignum-ref res len-x) (if minusp 1 0))))
972    (if fully-normalize
973      (%normalize-bignum-macro res)
974      (%mostly-normalize-bignum-macro res))))
975
976;;; NEGATE-BIGNUM-IN-PLACE -- Internal.
977;;;
978;;; This assumes bignum is positive; that is, the result of negating it will
979;;; stay in the provided allocated bignum.
980;;;
981(defun negate-bignum-in-place (bignum)
982  (bignum-negate-loop-really bignum (%bignum-length bignum) bignum)
983  bignum)
984
985
986 
987
988(defun copy-bignum (bignum)
989  (let ((res (%allocate-bignum (%bignum-length bignum))))
990    (bignum-replace res bignum)
991    res))
992
993
994
995;;; BIGNUM-ASHIFT-RIGHT -- Public.
996;;;
997;;; First compute the number of whole digits to shift, shifting them by
998;;; skipping them when we start to pick up bits, and the number of bits to
999;;; shift the remaining digits into place.  If the number of digits is greater
1000;;; than the length of the bignum, then the result is either 0 or -1.  If we
1001;;; shift on a digit boundary (that is, n-bits is zero), then we just copy
1002;;; digits.  The last branch handles the general case which uses a macro that a
1003;;; couple other routines use.  The fifth argument to the macro references
1004;;; locals established by the macro.
1005;;;
1006
1007
1008(defun bignum-ashift-right (bignum x)
1009  (declare (type bignum-type bignum)
1010           (fixnum x)
1011           (optimize (speed 3) (safety 0)))
1012  (let ((bignum-len (%bignum-length bignum)))
1013    (declare (type bignum-index bignum-len))
1014    (multiple-value-bind (digits n-bits) (truncate x digit-size)
1015      (declare (type bignum-index digits)(fixnum n-bits))
1016      (cond
1017       ((>= digits bignum-len)
1018        (if (bignum-plusp bignum) 0 -1))
1019       ((eql 0 n-bits)
1020        (bignum-ashift-right-digits bignum digits))
1021       (t
1022        (shift-right-unaligned bignum digits n-bits (- bignum-len digits)
1023                                      ((= j res-len-1)
1024                                       (setf (bignum-ref res j)
1025                                             (%ashr (bignum-ref bignum i) n-bits))
1026                                       (%normalize-bignum-macro res))
1027                                      res))))))
1028
1029                               
1030
1031
1032
1033;;; BIGNUM-ASHIFT-RIGHT-DIGITS -- Internal.
1034;;;
1035(defun bignum-ashift-right-digits (bignum digits)
1036  (declare (type bignum-type bignum)
1037           (type bignum-index digits))
1038  (let* ((res-len (- (%bignum-length bignum) digits))
1039         (res (%allocate-bignum res-len)))
1040    (declare (type bignum-index res-len)
1041             (type bignum-type res))
1042    (bignum-replace res bignum :start2 digits)
1043    (%normalize-bignum-macro res)))
1044
1045
1046;;; BIGNUM-BUFFER-ASHIFT-RIGHT -- Internal.
1047;;;
1048;;; GCD uses this for an in-place shifting operation.  This is different enough
1049;;; from BIGNUM-ASHIFT-RIGHT that it isn't worth folding the bodies into a
1050;;; macro, but they share the basic algorithm.  This routine foregoes a first
1051;;; test for digits being greater than or equal to bignum-len since that will
1052;;; never happen for its uses in GCD.  We did fold the last branch into a macro
1053;;; since it was duplicated a few times, and the fifth argument to it
1054;;; references locals established by the macro.
1055;;;
1056 
1057
1058;;; BIGNUM-ASHIFT-LEFT -- Public.
1059;;;
1060;;; This handles shifting a bignum buffer to provide fresh bignum data for some
1061;;; internal routines.  We know bignum is safe when called with bignum-len.
1062;;; First we compute the number of whole digits to shift, shifting them
1063;;; starting to store farther along the result bignum.  If we shift on a digit
1064;;; boundary (that is, n-bits is zero), then we just copy digits.  The last
1065;;; branch handles the general case.
1066;;;
1067(defun bignum-ashift-left (bignum x &optional bignum-len)
1068  (declare (type bignum-type bignum)
1069           (fixnum x)
1070           (type (or null bignum-index) bignum-len))
1071  (multiple-value-bind (digits n-bits)
1072                       (truncate x digit-size)
1073    (declare (fixnum digits n-bits))
1074    (let* ((bignum-len (or bignum-len (%bignum-length bignum)))
1075           (res-len (+ digits bignum-len 1)))
1076      (declare (fixnum bignum-len res-len))
1077      (when (> res-len maximum-bignum-length)
1078        (error "Can't represent result of left shift."))
1079      (if (zerop n-bits)
1080        (bignum-ashift-left-digits bignum bignum-len digits)
1081        (bignum-ashift-left-unaligned bignum digits n-bits res-len)))))
1082
1083;;; BIGNUM-ASHIFT-LEFT-DIGITS -- Internal.
1084;;;
1085(defun bignum-ashift-left-digits (bignum bignum-len digits)
1086  (declare (type bignum-index bignum-len digits))
1087  (let* ((res-len (+ bignum-len digits))
1088         (res (%allocate-bignum res-len)))
1089    (declare (type bignum-index res-len))
1090    (bignum-replace res bignum :start1 digits :end1 res-len :end2 bignum-len
1091                    :from-end t)
1092    res))
1093
1094
1095
1096;;; BIGNUM-ASHIFT-LEFT-UNALIGNED -- Internal.
1097;;;
1098;;; BIGNUM-TRUNCATE uses this to store into a bignum buffer by supplying res.
1099;;; When res comes in non-nil, then this foregoes allocating a result, and it
1100;;; normalizes the buffer instead of the would-be allocated result.
1101;;;
1102;;; We start storing into one digit higher than digits, storing a whole result
1103;;; digit from parts of two contiguous digits from bignum.  When the loop
1104;;; finishes, we store the remaining bits from bignum's first digit in the
1105;;; first non-zero result digit, digits.  We also grab some left over high
1106;;; bits from the last digit of bignum.
1107;;;
1108
1109(defun bignum-ashift-left-unaligned (bignum digits n-bits res-len
1110                                            &optional (res nil resp))
1111  (declare (type bignum-index digits res-len)
1112           (type (mod #.digit-size) n-bits))
1113  (let* ((remaining-bits (- digit-size n-bits))
1114         (res-len-1 (1- res-len))
1115         (res (or res (%allocate-bignum res-len))))
1116    (declare (type bignum-index res-len res-len-1)
1117             (optimize (speed 3) (safety 0)))
1118    (do ((i 0 i+1)
1119         (i+1 1 (1+ i+1))
1120         (j (1+ digits) (1+ j)))
1121        ((= j res-len-1)
1122         (setf (bignum-ref res digits)
1123               (%ashl (bignum-ref bignum 0) n-bits))
1124         (setf (bignum-ref res j)
1125               (%ashr (bignum-ref bignum i) remaining-bits))
1126         (if resp
1127           (%zero-trailing-sign-digits res res-len)
1128           (%mostly-normalize-bignum-macro res)))
1129      (declare (type bignum-index i i+1 j))
1130      (setf (bignum-ref res j)
1131            (%logior (%digit-logical-shift-right (bignum-ref bignum i)
1132                                                 remaining-bits)
1133                     (%ashl (bignum-ref bignum i+1) n-bits))))))
1134
1135
1136
1137
1138
1139
1140
1141;;;; Relational operators.
1142
1143
1144
1145;;; BIGNUM-COMPARE -- Public.
1146;;;
1147;;; This compares two bignums returning -1, 0, or 1, depending on whether a
1148;;; is less than, equal to, or greater than b.
1149;;;
1150;(proclaim '(function bignum-compare (bignum bignum) (integer -1 1)))
1151(defun bignum-compare (a b)
1152  (declare (type bignum-type a b))
1153  (let* ((a-plusp (bignum-plusp a))
1154         (b-plusp (bignum-plusp b)))
1155    (if (eq a-plusp b-plusp)
1156      (let* ((len-a (%bignum-length a))
1157             (len-b (%bignum-length b)))
1158        (declare (type bignum-index len-a len-b))
1159        (cond ((= len-a len-b)
1160               (do* ((i (1- len-a) (1- i)))
1161                    ((zerop i) (%compare-digits a b 0))
1162                 (declare (fixnum i))
1163                 (let* ((signum (%compare-digits a b i)))
1164                   (declare (fixnum signum))
1165                   (unless (zerop signum)
1166                     (return signum)))))
1167              ((> len-a len-b)
1168               (if a-plusp 1 -1))
1169              (t (if a-plusp -1 1))))
1170      (if a-plusp 1 -1))))
1171
1172
1173
1174
1175
1176
1177;;;; Integer length and logcount
1178
1179
1180(defun bignum-integer-length (big)
1181  (the fixnum (- (the fixnum (ash (the fixnum (%bignum-length big)) 5))
1182                 (the fixnum (%bignum-sign-bits big)))))
1183
1184; (not (zerop (logand integer1 integer2)
1185
1186(defun bignum-logtest (num1 num2)
1187  (let* ((length1 (%bignum-length num1))
1188         (length2 (%bignum-length num2))
1189         (n1-minusp (bignum-minusp num1))
1190         (n2-minusp (bignum-minusp num2)))
1191    (declare (fixnum length1 length2))
1192    (if (and n1-minusp n2-minusp) ; both neg, get out quick
1193      T       
1194      (or (dotimes (i (min length1 length2))
1195            (unless (zerop (the fixnum
1196                             (logand (the fixnum (bignum-ref num1 i))
1197                                     (the fixnum (bignum-ref num2 i)))))
1198              (return t)))
1199          (if (< length1 length2)
1200            n1-minusp
1201            (if (< length1 length2)
1202              n2-minusp))))))
1203
1204(defun logtest-fix-big (fix big)
1205  (declare (fixnum fix))
1206  (unless (zerop fix)
1207    (if (plusp fix)
1208      (or
1209       (not (eql 0 (the fixnum (logand (the fixnum (bignum-ref big 0)) fix))))
1210       (and (> (%bignum-length big) 1)
1211            (not (eql 0 (the fixnum (logand (the fixnum (bignum-ref big 1))
1212                                            (the fixnum (ash fix -32))))))))
1213      t)))
1214
1215
1216(defun bignum-logcount (bignum)
1217  (declare (type bignum-type bignum))
1218  (let* ((length (%bignum-length bignum))
1219         (plusp (bignum-plusp bignum))
1220         (result 0))
1221    (declare (type bignum-index length)
1222             (fixnum result))
1223    (if plusp
1224      (dotimes (index length result)
1225        (incf result (the fixnum (%logcount bignum index))))
1226      (dotimes (index length result)
1227        (incf result (the fixnum (%logcount-complement bignum index)))))))
1228
1229
1230;;;; Logical operations.
1231
1232;;; NOT.
1233;;;
1234
1235;;; BIGNUM-LOGICAL-NOT -- Public.
1236;;;
1237(defun bignum-logical-not (a)
1238  (declare (type bignum-type a))
1239  (let* ((len (%bignum-length a))
1240         (res (%allocate-bignum len)))
1241    (declare (type bignum-index len))
1242    (dotimes (i len res)
1243      (bignum-set res i (%lognot (the fixnum (bignum-ref a i)))))))
1244
1245
1246
1247
1248;;; AND.
1249;;;
1250
1251;;; BIGNUM-LOGICAL-AND -- Public.
1252;;;
1253(defun bignum-logical-and (a b)
1254  (declare (type bignum-type a b))
1255  (let* ((len-a (%bignum-length a))
1256         (len-b (%bignum-length b))
1257         (shorter a)
1258         (longer b)
1259         (shorter-len len-a)
1260         (longer-len len-b)
1261         (shorter-positive (bignum-plusp a)))
1262    (declare (type bignum-index len-a len-b shorter-len longer-len))
1263    (when (< len-b len-a)
1264      (setq shorter b
1265            longer a
1266            shorter-len len-b
1267            longer-len len-a
1268            shorter-positive (bignum-plusp b)))
1269    (let* ((result (%allocate-bignum longer-len)))
1270      (%bignum-logand shorter-len shorter longer result)
1271      (unless shorter-positive
1272        (bignum-replace result longer :start1 shorter-len :start2 shorter-len :end1 longer-len :end2 longer-len))
1273      (%normalize-bignum-macro result))))
1274
1275
1276;;;
1277;;;
1278;;; bignum-logandc2
1279
1280(defun bignum-logandc2 (a b)
1281  (declare (type bignum-type a b))
1282  (let* ((len-a (%bignum-length a))
1283         (len-b (%bignum-length b))
1284         (a-plusp (bignum-plusp a))
1285         (b-plusp (bignum-plusp b)))
1286    (declare (type bignum-index len-a len-b))
1287    (cond
1288     ((< len-a len-b)
1289      (logandc2-shorter-any a len-a b len-b (if a-plusp (%allocate-bignum len-a) (%allocate-bignum len-b))))
1290     ((< len-b len-a) ; b shorter
1291      (logandc1-shorter-any b len-b a len-a (if b-plusp (%allocate-bignum len-a)(%allocate-bignum len-b))))
1292     (t (logandc2-shorter-any a len-a b len-b (%allocate-bignum len-a))))))
1293
1294(defun logandc2-shorter-any (a len-a b len-b res)
1295  (declare (type bignum-type a b res)
1296           (type bignum-index len-a len-b))
1297  (dotimes (i len-a)
1298    (setf (bignum-ref res i)
1299          (logand (the fixnum (bignum-ref a i))
1300                  (the fixnum (%lognot (the fixnum (bignum-ref b i)))))))
1301  (if (bignum-minusp a)
1302    (do ((i len-a (1+ i)))
1303          ((= i len-b))
1304        (declare (type bignum-index i))
1305      (setf (bignum-ref res i)
1306            (%lognot (the fixnum (bignum-ref b i))))))
1307  (%normalize-bignum-macro res))
1308
1309
1310
1311(defun logandc1-shorter-any (a len-a b len-b res)
1312  (declare (type bignum-type a b res)
1313           (type bignum-index len-a len-b))
1314  (dotimes (i len-a)
1315    (setf (bignum-ref res i)
1316          (logand
1317           (the fixnum (%lognot (the fixnum (bignum-ref a i))))
1318           (the fixnum (bignum-ref b i)))))
1319  (when (bignum-plusp a)
1320    (unless (= len-a len-b)
1321      (bignum-replace res b :start1 len-a :start2 len-a :end1 len-b :end2 len-b)))
1322  (%normalize-bignum-macro res))
1323
1324
1325
1326(defun fix-big-logand (fix big)
1327  (let* ((len-b (%bignum-length big))
1328         (res (if (< fix 0)(%allocate-bignum len-b))))
1329    (declare (fixnum fix len-b))       
1330    (let ((val (fix-digit-logand fix big res)))
1331      (if res
1332        (progn
1333          (bignum-replace res big :start1 2 :start2 2 :end1 len-b :end2 len-b)
1334          (%normalize-bignum-macro res))
1335        val))))
1336
1337
1338(defun fix-big-logandc2 (fix big)
1339  (let* ((len-b (%bignum-length big))
1340         (res (if (< fix 0)(%allocate-bignum len-b))))
1341    (declare (fixnum fix len-b))       
1342    (let ((val (fix-digit-logandc2 fix big res)))
1343      (if res
1344        (progn
1345          (do ((i 2 (1+ i)))
1346              ((= i len-b))
1347            (declare (type bignum-index i))
1348            (setf (bignum-ref res i)
1349                  (%lognot (bignum-ref big i))))
1350          (%normalize-bignum-macro res))
1351        val))))
1352
1353(defun fix-big-logandc1 (fix big)
1354  (let* ((len-b (%bignum-length big))
1355         (res (if (>= fix 0)(%allocate-bignum len-b))))
1356    (declare (fixnum fix len-b))       
1357    (let ((val (fix-digit-logandc1 fix big res)))
1358      (if res
1359        (progn 
1360          (bignum-replace res big :start1 2 :start2 2 :end1 len-b :end2 len-b)
1361          (%normalize-bignum-macro res))
1362        val))))
1363
1364
1365;;; IOR.
1366;;;
1367
1368;;; BIGNUM-LOGICAL-IOR -- Public.
1369;;;
1370(defun bignum-logical-ior (a b)
1371  (declare (type bignum-type a b))
1372  (let* ((len-a (%bignum-length a))
1373         (len-b (%bignum-length b))
1374         (longer-len len-b)
1375         (shorter-len len-a)
1376         (shorter a)
1377         (longer b)
1378         (shorter-positive (bignum-plusp a)))
1379    (declare (type bignum-index len-a len-b longer-len shorter-len))
1380    (when (< len-b len-a)
1381      (setq shorter b
1382            longer a
1383            shorter-len len-b
1384            longer-len len-a
1385            shorter-positive (bignum-plusp b)))
1386    (let* ((result (%allocate-bignum longer-len)))
1387      (%bignum-logior shorter-len shorter longer result)
1388      (unless (= shorter-len longer-len)
1389        (if shorter-positive
1390          (bignum-replace result longer :start1 shorter-len :start2 shorter-len :end1 longer-len :end2 longer-len)
1391          (do* ((i shorter-len (1+ i)))
1392               ((= i longer-len))
1393            (declare (type bignum-index i))
1394            (setf (bignum-ref result i) #xffffffff))))
1395      (%normalize-bignum-macro result))))
1396
1397
1398
1399;;; XOR.
1400;;;
1401
1402;;; BIGNUM-LOGICAL-XOR -- Public.
1403;;;
1404(defun bignum-logical-xor (a b)
1405  (declare (type bignum-type a b))
1406  (let ((len-a (%bignum-length a))
1407        (len-b (%bignum-length b)))
1408    (declare (type bignum-index len-a len-b))
1409    (if (< len-a len-b)
1410        (bignum-logical-xor-aux a len-a b len-b (%allocate-bignum len-b))
1411        (bignum-logical-xor-aux b len-b a len-a (%allocate-bignum len-a)))))
1412
1413;;; BIGNUM-LOGICAL-XOR-AUX -- Internal.
1414;;;
1415;;; This takes the the shorter of two bignums in a and len-a.  Res is len-b
1416;;; long.  Do the XOR.
1417;;;
1418(defun bignum-logical-xor-aux (a len-a b len-b res)
1419  (declare (type bignum-type a b res)
1420           (type bignum-index len-a len-b))
1421  (dotimes (i len-a)
1422    (setf (bignum-ref res i)
1423          (%logxor (the fixnum (bignum-ref a i))
1424                  (the fixnum (bignum-ref b i)))))
1425  (unless (= len-a len-b)
1426    (let ((sign (if (bignum-minusp a) all-ones-digit 0)))
1427      (declare (fixnum sign))
1428      (do ((i len-a (1+ i)))
1429          ((= i len-b))
1430        (declare (type bignum-index i))
1431        (setf (bignum-ref res i)
1432              (%logxor (bignum-ref b i) sign)))))
1433  (%normalize-bignum-macro res))
1434
1435
1436;;;; TRUNCATE
1437
1438;;; Divide X by Y when Y is a single bignum digit. BIGNUM-TRUNCATE
1439;;; fixes up the quotient and remainder with respect to sign and
1440;;; normalization.
1441;;;
1442;;; We don't have to worry about shifting Y to make its most
1443;;; significant digit sufficiently large for %FLOOR to return
1444;;; digit-size quantities for the q-digit and r-digit. If Y is
1445;;; a single digit bignum, it is already large enough for
1446;;; %FLOOR. That is, it has some bits on pretty high in the
1447;;; digit.
1448
1449(defun bignum-truncate-single-digit (x len-x y)
1450  (declare (type bignum-index len-x))
1451  (let ((q (%allocate-bignum len-x))
1452        (r 0)
1453        (y (bignum-ref y 0)))
1454    (declare (type bignum-element-type r y))
1455    (do ((i (1- len-x) (1- i)))
1456        ((minusp i))
1457      (multiple-value-bind (q-digit r-digit)
1458          (%floor r (bignum-ref x i) y)
1459        (declare (type bignum-element-type q-digit r-digit))
1460        (setf (bignum-ref q i) q-digit)
1461        (setf r r-digit)))
1462    (let ((rem (%allocate-bignum 1)))
1463      (setf (bignum-ref rem 0) r)
1464      (values q rem))))
1465
1466;;; This returns a guess for the next division step. Y1 is the
1467;;; highest y digit, and y2 is the second to highest y
1468;;; digit. The x... variables are the three highest x digits
1469;;; for the next division step.
1470;;;
1471;;; From Knuth, our guess is either all ones or x-i and x-i-1
1472;;; divided by y1, depending on whether x-i and y1 are the
1473;;; same. We test this guess by determining whether guess*y2
1474;;; is greater than the three high digits of x minus guess*y1
1475;;; shifted left one digit:
1476;;;    ------------------------------
1477;;;   |    x-i    |   x-i-1  | x-i-2 |
1478;;;    ------------------------------
1479;;;    ------------------------------
1480;;; - | g*y1 high | g*y1 low |   0   |
1481;;;    ------------------------------
1482;;;             ...               <   guess*y2     ???   
1483;;; If guess*y2 is greater, then we decrement our guess by one
1484;;; and try again.  This returns a guess that is either
1485;;; correct or one too large.
1486(defun bignum-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
1487  (declare (type bignum-element-type y1 y2 x-i x-i-1 x-i-2))
1488  (let ((guess (if (= x-i y1)
1489                 all-ones-digit
1490                 (%floor x-i x-i-1 y1))))
1491    (declare (type bignum-element-type guess))
1492    (loop
1493      (multiple-value-bind (high-guess*y1 low-guess*y1)
1494          (%multiply guess y1)
1495        (declare (type bignum-element-type low-guess*y1
1496                       high-guess*y1))
1497        (multiple-value-bind (high-guess*y2 low-guess*y2)
1498            (%multiply guess y2)
1499          (declare (type bignum-element-type high-guess*y2
1500                         low-guess*y2))
1501          (multiple-value-bind (middle-digit borrow)
1502              (%subtract-with-borrow x-i-1 low-guess*y1 1)
1503            (declare (type bignum-element-type middle-digit)
1504                     (fixnum borrow))
1505            ;; Supplying borrow of 1 means there was no
1506            ;; borrow, and we know x-i-2 minus 0 requires
1507            ;; no borrow.
1508            (let ((high-digit (%subtract-with-borrow x-i
1509                                                     high-guess*y1
1510                                                     borrow)))
1511              (declare (type bignum-element-type high-digit))
1512              (if (and (= high-digit 0)
1513                       (or (> high-guess*y2
1514                              middle-digit)
1515                           (and (= middle-digit
1516                                   high-guess*y2)
1517                                (> low-guess*y2
1518                                   x-i-2))))
1519                (setf guess (%subtract-with-borrow guess 1 1))
1520                (return guess)))))))))
1521
1522
1523;;; This returns the amount to shift y to place a one in the
1524;;; second highest bit. Y must be positive. If the last digit
1525;;; of y is zero, then y has a one in the previous digit's
1526;;; sign bit, so we know it will take one less than digit-size
1527;;; to get a one where we want. Otherwise, we count how many
1528;;; right shifts it takes to get zero; subtracting this value
1529;;; from digit-size tells us how many high zeros there are
1530;;; which is one more than the shift amount sought.
1531;;;
1532;;; Note: This is exactly the same as one less than the
1533;;; integer-length of the last digit subtracted from the
1534;;; digit-size.
1535;;;
1536;;; We shift y to make it sufficiently large that doing the
1537;;; 2*digit-size by digit-size %FLOOR calls ensures the quotient and
1538;;; remainder fit in digit-size.
1539(defun shift-y-for-truncate (y)
1540  (the fixnum (1- (the fixnum (%bignum-sign-bits y)))))
1541
1542;;; Stores two bignums into the truncation bignum buffers,
1543;;; shifting them on the way in. This assumes x and y are
1544;;; positive and at least two in length, and it assumes
1545;;; truncate-x and truncate-y are one digit longer than x and
1546;;; y.
1547(defun shift-and-store-truncate-buffers (truncate-x truncate-y x len-x y len-y shift)
1548  (declare (type bignum-index len-x len-y)
1549           (type (integer 0 (#.digit-size)) shift))
1550  (cond ((zerop shift)
1551         (bignum-replace truncate-x x :end1 len-x)
1552         (bignum-replace truncate-y y :end1 len-y))
1553        (t
1554         (bignum-ashift-left-unaligned x 0 shift (1+ len-x)
1555                                       truncate-x)
1556         (bignum-ashift-left-unaligned y 0 shift (1+ len-y)
1557                                       truncate-y))))
1558
1559;;; Divide TRUNCATE-X by TRUNCATE-Y, returning the quotient
1560;;; and destructively modifying TRUNCATE-X so that it holds
1561;;; the remainder.
1562;;;
1563;;; LEN-X and LEN-Y tell us how much of the buffers we care about.
1564;;;
1565;;; TRUNCATE-X definitely has at least three digits, and it has one
1566;;; more than TRUNCATE-Y. This keeps i, i-1, i-2, and low-x-digit
1567;;; happy. Thanks to SHIFT-AND-STORE-TRUNCATE-BUFFERS.
1568
1569(defun do-truncate (truncate-x truncate-y len-x len-y)
1570  (declare (type bignum-index len-x len-y))
1571  (let* ((len-q (- len-x len-y))
1572         ;; Add one for extra sign digit in case high bit is on.
1573         (q (%allocate-bignum (1+ len-q)))
1574         (k (1- len-q))
1575         (y1 (bignum-ref truncate-y (1- len-y)))
1576         (y2 (bignum-ref truncate-y (- len-y 2)))
1577         (i (1- len-x))
1578         (i-1 (1- i))
1579         (i-2 (1- i-1))
1580         (low-x-digit (- i len-y)))
1581    (declare (type bignum-index len-q k i i-1 i-2 low-x-digit)
1582             (type bignum-element-type y1 y2))
1583    (loop
1584      (setf (bignum-ref q k)
1585            (try-bignum-truncate-guess
1586             truncate-x truncate-y
1587             ;; This modifies TRUNCATE-X. Must access
1588             ;; elements each pass.
1589             (bignum-truncate-guess y1 y2
1590                                    (bignum-ref truncate-x i)
1591                                    (bignum-ref truncate-x i-1)
1592                                    (bignum-ref truncate-x i-2))
1593             len-y low-x-digit))
1594      (cond ((zerop k) (return))
1595            (t (decf k)
1596               (decf low-x-digit)
1597               (shiftf i i-1 i-2 (1- i-2)))))
1598    q))
1599
1600#+notyet
1601(defun do-truncate-no-quo (truncate-x truncate-y len-x len-y)
1602  (declare (type bignum-index len-x len-y))
1603  (let* ((len-q (- len-x len-y))
1604         (k (1- len-q))
1605         (i (1- len-x))
1606         (low-x-digit (- i len-y)))
1607    (declare (type bignum-index len-q k i  low-x-digit))
1608    (loop
1609      (let* ((guess (bignum-truncate-guess truncate-x i truncate-y (the fixnum (1- len-y)))                                 
1610        (try-bignum-truncate-guess guess len-y low-x-digit)
1611        (cond ((zerop k) (return))
1612              (t (decf k)
1613                 (decf low-x-digit)
1614                 (setq i (1- i))))))
1615    nil))))
1616
1617;;; This takes a digit guess, multiplies it by TRUNCATE-Y for a
1618;;; result one greater in length than LEN-Y, and subtracts this result
1619;;; from TRUNCATE-X. LOW-X-DIGIT is the first digit of X to start
1620;;; the subtraction, and we know X is long enough to subtract a LEN-Y
1621;;; plus one length bignum from it. Next we check the result of the
1622;;; subtraction, and if the high digit in X became negative, then our
1623;;; guess was one too big. In this case, return one less than GUESS
1624;;; passed in, and add one value of Y back into X to account for
1625;;; subtracting one too many. Knuth shows that the guess is wrong on
1626;;; the order of 3/b, where b is the base (2 to the digit-size power)
1627;;; -- pretty rarely.
1628
1629(defun try-bignum-truncate-guess (truncate-x truncate-y guess len-y low-x-digit)
1630  (declare (type bignum-index low-x-digit len-y)
1631           (type bignum-element-type guess))
1632  (let ((carry-digit 0)
1633        (borrow 1)
1634        (i low-x-digit))
1635    (declare (type bignum-element-type carry-digit)
1636             (type bignum-index i)
1637             (fixnum borrow))
1638    ;; Multiply guess and divisor, subtracting from dividend
1639    ;; simultaneously.
1640    (dotimes (j len-y)
1641      (multiple-value-bind (high-digit low-digit)
1642          (%multiply-and-add3 guess
1643                              (bignum-ref truncate-y j)
1644                              carry-digit)
1645        (declare (type bignum-element-type high-digit low-digit))
1646        (setf carry-digit high-digit)
1647        (multiple-value-bind (x temp-borrow)
1648            (%subtract-with-borrow (bignum-ref truncate-x i)
1649                                   low-digit
1650                                   borrow)
1651          (declare (type bignum-element-type x)
1652                   (fixnum temp-borrow))
1653          (setf (bignum-ref truncate-x i) x)
1654          (setf borrow temp-borrow)))
1655      (incf i))
1656    (setf (bignum-ref truncate-x i)
1657          (%subtract-with-borrow (bignum-ref truncate-x i)
1658                                 carry-digit borrow))
1659    ;; See whether guess is off by one, adding one
1660    ;; Y back in if necessary.
1661    (cond ((%digit-0-or-plusp (bignum-ref truncate-x i))
1662           guess)
1663          (t
1664           ;; If subtraction has negative result, add one
1665           ;; divisor value back in. The guess was one too
1666           ;; large in magnitude.
1667           (let ((i low-x-digit)
1668                 (carry 0))
1669             (dotimes (j len-y)
1670               (multiple-value-bind (v k)
1671                   (%add-with-carry (bignum-ref truncate-y j)
1672                                    (bignum-ref truncate-x i)
1673                                    carry)
1674                 (declare (type bignum-element-type v))
1675                 (setf (bignum-ref truncate-x i) v)
1676                 (setf carry k))
1677               (incf i))
1678             (setf (bignum-ref truncate-x i)
1679                   (%add-with-carry (bignum-ref truncate-x i)
1680                                    0 carry)))
1681           (%subtract-with-borrow guess 1 1)))))
1682
1683;;; Someone (from the original CMUCL or SPICE Lisp project, perhaps)
1684;;; is the "I" who implemented the original version of this.
1685
1686;;; This is the original sketch of the algorithm from which I implemented this
1687;;; TRUNCATE, assuming both operands are bignums. I should modify this to work
1688;;; with the documentation on my functions, as a general introduction. I've
1689;;; left this here just in case someone needs it in the future. Don't look at
1690;;; this unless reading the functions' comments leaves you at a loss. Remember
1691;;; this comes from Knuth, so the book might give you the right general
1692;;; overview.
1693;;;
1694;;; (truncate x y):
1695;;;
1696;;; If X's magnitude is less than Y's, then result is 0 with remainder X.
1697;;;
1698;;; Make x and y positive, copying x if it is already positive.
1699;;;
1700;;; Shift y left until there's a 1 in the 30'th bit (most significant, non-sign
1701;;;       digit)
1702;;;    Just do most sig digit to determine how much to shift whole number.
1703;;; Shift x this much too.
1704;;; Remember this initial shift count.
1705;;;
1706;;; Allocate q to be len-x minus len-y quantity plus 1.
1707;;;
1708;;; i = last digit of x.
1709;;; k = last digit of q.
1710;;;
1711;;; LOOP
1712;;;
1713;;; j = last digit of y.
1714;;;
1715;;; compute guess.
1716;;; if x[i] = y[j] then g = (1- (ash 1 digit-size))
1717;;; else g = x[i]x[i-1]/y[j].
1718;;;
1719;;; check guess.
1720;;; %UNSIGNED-MULTIPLY returns b and c defined below.
1721;;;    a = x[i-1] - (logand (* g y[j]) #xFFFFFFFF).
1722;;;       Use %UNSIGNED-MULTIPLY taking low-order result.
1723;;;    b = (logand (ash (* g y[j-1]) (- digit-size)) (1- (ash 1 digit-size))).
1724;;;    c = (logand (* g y[j-1]) (1- (ash 1 digit-size))).
1725;;; if a < b, okay.
1726;;; if a > b, guess is too high
1727;;;    g = g - 1; go back to "check guess".
1728;;; if a = b and c > x[i-2], guess is too high
1729;;;    g = g - 1; go back to "check guess".
1730;;; GUESS IS 32-BIT NUMBER, SO USE THING TO KEEP IN SPECIAL REGISTER
1731;;; SAME FOR A, B, AND C.
1732;;;
1733;;; Subtract g * y from x[i - len-y+1]..x[i]. See paper for doing this in step.
1734;;; If x[i] < 0, guess is screwed up.
1735;;;    negative g, then add 1
1736;;;    zero or positive g, then subtract 1
1737;;; AND add y back into x[len-y+1..i].
1738;;;
1739;;; q[k] = g.
1740;;; i = i - 1.
1741;;; k = k - 1.
1742;;;
1743;;; If k>=0, goto LOOP.
1744;;;
1745;;; Now quotient is good, but remainder is not.
1746;;; Shift x right by saved initial left shifting count.
1747;;;
1748;;; Check quotient and remainder signs.
1749;;; x pos y pos --> q pos r pos
1750;;; x pos y neg --> q neg r pos
1751;;; x neg y pos --> q neg r neg
1752;;; x neg y neg --> q pos r neg
1753;;;
1754;;; Normalize quotient and remainder. Cons result if necessary.
1755
1756
1757(defun bignum-truncate (x y &optional no-rem)
1758  (declare (type bignum-type x y))
1759  (DECLARE (IGNORE NO-REM))
1760  ;; Divide X by Y returning the quotient and remainder. In the
1761  ;; general case, we shift Y to set up for the algorithm, and we
1762  ;; use two buffers to save consing intermediate values. X gets
1763  ;; destructively modified to become the remainder, and we have
1764  ;; to shift it to account for the initial Y shift. After we
1765  ;; multiple bind q and r, we first fix up the signs and then
1766  ;; return the normalized results.
1767  (let* ((x-plusp (%bignum-0-or-plusp x (%bignum-length x)))
1768         (y-plusp (%bignum-0-or-plusp y (%bignum-length y)))
1769         (x (if x-plusp x (negate-bignum x nil)))
1770         (y (if y-plusp y (negate-bignum y nil)))
1771         (len-x (%bignum-length x))
1772         (len-y (%bignum-length y)))
1773    (multiple-value-bind (q r)
1774        (cond ((< len-y 2)
1775               (bignum-truncate-single-digit x len-x y))
1776              ((plusp (bignum-compare y x))
1777               (let ((res (%allocate-bignum len-x)))
1778                 (dotimes (i len-x)
1779                   (setf (bignum-ref res i) (bignum-ref x i)))
1780                 (values 0 res)))
1781              (t
1782               (let ((len-x+1 (1+ len-x)))
1783                 (with-bignum-buffers ((truncate-x len-x+1)
1784                                       (truncate-y (1+ len-y)))
1785                   (let ((y-shift (shift-y-for-truncate y)))
1786                     (shift-and-store-truncate-buffers truncate-x
1787                                                       truncate-y
1788                                                       x len-x
1789                                                       y len-y
1790                                                       y-shift)
1791                     (values
1792                      (do-truncate truncate-x
1793                        truncate-y
1794                        len-x+1
1795                        len-y)
1796                      ;; Now DO-TRUNCATE has executed, we just
1797                      ;; tidy up the remainder (in TRUNCATE-X)
1798                      ;; and return it.
1799                      (cond
1800                        ((zerop y-shift)
1801                         (let ((res (%allocate-bignum len-y)))
1802                           (declare (type bignum-type res))
1803                           (bignum-replace res truncate-x :end2 len-y)
1804                           (%normalize-bignum-macro res)))
1805                        (t
1806                         (shift-right-unaligned
1807                          truncate-x 0 y-shift len-y
1808                          ((= j res-len-1)
1809                           (setf (bignum-ref res j)
1810                                 (%ashr (bignum-ref truncate-x i)
1811                                        y-shift))
1812                           (%normalize-bignum-macro res))
1813                          res)))))))))
1814      (let ((quotient (cond ((eq x-plusp y-plusp) q)
1815                            ((typep q 'fixnum) (the fixnum (- q)))
1816                            (t (negate-bignum-in-place q))))
1817            (rem (cond (x-plusp r)
1818                       ((typep r 'fixnum) (the fixnum (- r)))
1819                       (t (negate-bignum-in-place r)))))
1820        (values (if (typep quotient 'fixnum)
1821                  quotient
1822                  (%normalize-bignum-macro quotient))
1823                (if (typep rem 'fixnum)
1824                  rem
1825                  (%normalize-bignum-macro rem)))))))
1826
1827(defun bignum-truncate-by-fixnum (bignum fixnum)
1828  (with-small-bignum-buffers ((y fixnum))
1829    (bignum-truncate bignum y)))
1830
1831(defun bignum-truncate-by-fixnum-no-quo (bignum fixnum)
1832  (nth-value 1 (bignum-truncate-by-fixnum bignum fixnum)))
1833
1834;;; This may do unnecessary computation in some cases.
1835(defun bignum-rem (x y)
1836  (nth-value 1 (bignum-truncate x y)))
1837
1838
1839
1840;;;; General utilities.
1841
1842(defun %zero-trailing-sign-digits (bignum len)
1843  (declare (fixnum len))
1844  (unless (<= len 1)
1845    (do ((next (bignum-ref bignum (the fixnum (- len 2)))
1846               (bignum-ref bignum (the fixnum (- len 2))))
1847         (sign (bignum-ref bignum (the fixnum (- len 1)))
1848               next))
1849        ((not (zerop (the fixnum (%logxor sign (%ashr next 31))))))
1850      (decf len)
1851      (setf (bignum-ref bignum len) 0)
1852      ;; Return, unless we've already done so (having found significant
1853      ;; digits earlier.)
1854      (when (= len 1)
1855        (return))))
1856  len)
1857
1858
1859(defun %normalize-bignum-2 (return-fixnum-p bignum)
1860  (let* ((len (%bignum-length bignum))
1861         (newlen (%zero-trailing-sign-digits bignum len)))
1862    (declare (fixnum len newlen))
1863    (unless (= len newlen)
1864      (%set-bignum-length newlen bignum))
1865    (or (and return-fixnum-p
1866             (%maybe-fixnum-from-one-or-two-digit-bignum bignum))
1867        bignum)))
1868           
1869   
1870;;; %MOSTLY-NORMALIZE-BIGNUM -- Internal.
1871;;;
1872;;; This drops the last digit if it is unnecessary sign information.  It
1873;;; repeats this as needed, possibly ending with a fixnum magnitude but never
1874;;; returning a fixnum.
1875;;;
1876
1877(defun %mostly-normalize-bignum (res &optional len)
1878  (declare (ignore len))
1879  (%normalize-bignum-2 nil res))
1880
1881
1882
1883
1884
1885(defun load-byte (size position integer)
1886  (if (and (bignump integer)
1887           (<= size (- 63 target::fixnumshift))
1888           (fixnump position))
1889    (%ldb-fixnum-from-bignum integer size position)
1890    (let ((mask (byte-mask size)))
1891      (if (and (fixnump mask) (fixnump integer)(fixnump position))
1892        (%ilogand mask (%iasr position integer))
1893        (logand mask (ash integer (- position)))))))
1894
1895
1896#+safe-but-slow
1897;;; This is basically the same algorithm as the "destructive"
1898;;; version below; while it may be more readable, it's often
1899;;; slower and conses too much to be at all viable.
1900(defun %bignum-bignum-gcd (u v)
1901  (setq u (abs u) v (abs v))
1902  (do* ((g 1 (ash g 1)))
1903       ((or (oddp u) (oddp v))
1904        (do* ()
1905             ((zerop u) (* g v))
1906          (cond ((evenp u) (setq u (ash u -1)))
1907                ((evenp v) (setq v (ash v -1)))
1908                (t (let* ((temp (ash (abs (- u v)) -1)))
1909                     (if (< u v)
1910                       (setq v temp)
1911                       (setq u temp)))))))
1912    (setq u (ash u -1) v (ash v -1))))
1913
1914
1915
1916
1917#-safe-but-slow
1918(progn
1919(defun %positive-bignum-bignum-gcd (u0 v0)
1920  (let* ((u-len (%bignum-length u0))
1921         (v-len (%bignum-length v0)))
1922    (declare (fixnum u-len v-len))
1923    (if (or (< u-len v-len)
1924            (and (= u-len v-len)
1925                 (< (bignum-compare u0 v0) 0)))
1926      (psetq u0 v0 v0 u0 u-len v-len v-len u-len))
1927    (with-bignum-buffers ((u u-len)
1928                          (u2 u-len)
1929                          (v v-len)
1930                          (v2 v-len))
1931      (bignum-replace u u0)
1932      (bignum-replace v v0)
1933      (let* ((u-trailing-0-bits (%bignum-count-trailing-zero-bits u))
1934             (u-trailing-0-digits (ash u-trailing-0-bits -5))
1935             (v-trailing-0-bits (%bignum-count-trailing-zero-bits v))
1936             (v-trailing-0-digits (ash v-trailing-0-bits -5)))
1937        (declare (fixnum u-trailing-0-bits v-trailing-0-bits))
1938        (unless (zerop u-trailing-0-bits)
1939          (bignum-shift-right-loop-1
1940           (logand u-trailing-0-bits 31)
1941           u2
1942           u
1943           (the fixnum (1- (the fixnum (- u-len u-trailing-0-digits ))))
1944           u-trailing-0-digits)
1945          (rotatef u u2)
1946          (%mostly-normalize-bignum-macro u)
1947          (setq u-len (%bignum-length u)))
1948        (unless (zerop v-trailing-0-bits)
1949          (bignum-shift-right-loop-1
1950           (logand v-trailing-0-bits 31)
1951           v2
1952           v
1953           (the fixnum (1- (the fixnum (- v-len v-trailing-0-digits))))
1954           v-trailing-0-digits)
1955          (rotatef v v2)
1956          (%mostly-normalize-bignum-macro v)
1957          (setq v-len (%bignum-length v)))
1958        (let* ((shift (min u-trailing-0-bits
1959                           v-trailing-0-bits)))
1960          (loop
1961              (let* ((fix-u (and (<= u-len 2)
1962                                 (%maybe-fixnum-from-one-or-two-digit-bignum u)))
1963                     (fix-v (and (<= v-len 2)
1964                                 (%maybe-fixnum-from-one-or-two-digit-bignum v))))
1965                (if fix-v
1966                  (if fix-u
1967                    (return (ash (%fixnum-gcd fix-u fix-v) shift))
1968                    (return (ash (bignum-fixnum-gcd u fix-v) shift)))
1969                  (if fix-u
1970                    (return (ash (bignum-fixnum-gcd v fix-u) shift)))))
1971              (let* ((signum (if (> u-len v-len)
1972                               1
1973                               (if (< u-len v-len)
1974                                 -1
1975                                 (bignum-compare u v)))))
1976                (declare (fixnum signum))
1977                (case signum
1978                  (0                    ; (= u v)
1979                   (if (zerop shift)
1980                     (let* ((copy (%allocate-bignum u-len)))
1981                       (bignum-replace copy u)
1982                       (return copy))
1983                     (return (ash u shift))))
1984                  (1                    ; (> u v)
1985                   (bignum-subtract-loop u u-len v v-len u)
1986                   (%mostly-normalize-bignum-macro u)
1987                   (setq u-len (%bignum-length u))
1988                   (setq u-trailing-0-bits
1989                         (%bignum-count-trailing-zero-bits u)
1990                         u-trailing-0-digits
1991                         (ash u-trailing-0-bits -5))
1992                   (unless (zerop u-trailing-0-bits)
1993                     (%init-misc 0 u2)
1994                     (bignum-shift-right-loop-1
1995                      (logand u-trailing-0-bits 31)
1996                      u2
1997                      u
1998                      (the fixnum (1- (the fixnum (- u-len
1999                                                     u-trailing-0-digits))))
2000                      u-trailing-0-digits)
2001                     (rotatef u u2)
2002                     (%mostly-normalize-bignum-macro u)
2003                     (setq u-len (%bignum-length u))))
2004                  (t                    ; (> v u)
2005                   (bignum-subtract-loop v v-len u u-len v)
2006                   (%mostly-normalize-bignum-macro v)
2007                   (setq v-len (%bignum-length v))
2008                   (setq v-trailing-0-bits
2009                         (%bignum-count-trailing-zero-bits v)
2010                         v-trailing-0-digits
2011                         (ash v-trailing-0-bits -5))
2012                   (unless (zerop v-trailing-0-bits)
2013                     (%init-misc 0 v2)
2014                     (bignum-shift-right-loop-1
2015                      (logand v-trailing-0-bits 31)
2016                      v2
2017                      v
2018                      (the fixnum (1- (the fixnum (- v-len v-trailing-0-digits))))
2019                      v-trailing-0-digits)
2020                     (rotatef v v2)
2021                     (%mostly-normalize-bignum-macro v)
2022                     (setq v-len (%bignum-length v))))))))))))
2023
2024(defun %bignum-bignum-gcd (u v)
2025  (with-negated-bignum-buffers u v %positive-bignum-bignum-gcd))
2026)
2027
2028
2029(defun bignum-shift-right-loop-1 (nbits result source len idx)
2030  (declare (type bignum-type result source)
2031           (type (mod 32) nbits)
2032           (type bignum-index idx len))
2033  (let* ((rbits (- 32 nbits)))
2034    (declare (type (mod 33) rbits))
2035    (dotimes (j len)
2036      (let* ((x (bignum-ref source idx)))
2037        (declare (type bignum-element-type x))
2038        (setq x (%ilsr nbits x))
2039        (incf idx)
2040        (let* ((y (bignum-ref source idx)))
2041          (declare (type bignum-element-type y))
2042          (setq y (%ashl y rbits))
2043          (setf (bignum-ref result j)
2044                (%logior x y)))))
2045    (setf (bignum-ref result len)
2046          (%ilsr nbits (bignum-ref source idx)))
2047    idx))
2048   
2049
2050(defun %logcount (bignum idx)
2051  (%ilogcount (bignum-ref bignum idx)))
2052
2053(defun %logcount-complement (bignum idx)
2054  (- 32 (the fixnum (%ilogcount (bignum-ref bignum idx)))))
2055
2056(defun %bignum-evenp (bignum)
2057  (not (logbitp 0 (the fixnum (bignum-ref bignum 0)))))
2058
2059(defun %bignum-oddp (bignum)
2060  (logbitp 0 (the fixnum (bignum-ref bignum 0))))
2061
2062(defun %ldb-fixnum-from-bignum (bignum size position)
2063  (declare (fixnum size position))
2064  (let* ((low-idx (ash position -5))
2065         (low-bit (logand position 31))
2066         (hi-bit (+ low-bit size))
2067         (len (%bignum-length bignum))
2068         (minusp (bignum-minusp bignum)))
2069    (declare (fixnum size position low-bit hi-bit low-idx len))
2070    (if (>= low-idx len)
2071      (if minusp (1- (ash 1 size)) 0)
2072      (flet ((ldb32 (digit size pos)
2073               (declare (fixnum digit size pos))
2074               (logand (the fixnum (1- (ash 1 size)))
2075                       (the fixnum (ash digit (the fixnum (- pos)))))))
2076        (let* ((low-digit (bignum-ref bignum low-idx))
2077               (chunk-lo (ldb32 low-digit (min size (%i- 32 low-bit)) low-bit)))
2078          (if (< hi-bit 32) 
2079            chunk-lo
2080            (let* ((have (- 32 low-bit))
2081                   (remain (- size have)))
2082              (declare (fixnum have remain))
2083              (setq low-idx (1+ low-idx))
2084              (when (> remain 32)
2085                (setq chunk-lo
2086                      (logior (ash (if (< low-idx len)
2087                                     (bignum-ref bignum low-idx)
2088                                     (if minusp all-ones-digit 0))
2089                                   have)
2090                              chunk-lo))
2091                (incf have 32)
2092                (decf remain 32)
2093                (incf low-idx))
2094              (let* ((high-digit
2095                      (if (>= low-idx len)
2096                        (if minusp all-ones-digit 0)
2097                        (bignum-ref bignum low-idx)))
2098                     (chunk-hi (ldb32 high-digit remain 0)))
2099                (%ilogior (ash chunk-hi have) chunk-lo)))))))))
2100
2101
2102
2103(defun bignum-negate-loop-really (big len res)
2104  (declare (fixnum len))
2105  (let* ((carry 1))
2106    (dotimes (i len carry)
2107      (multiple-value-bind (result-digit carry-out)
2108          (%add-with-carry (%lognot (bignum-ref big i)) 0 carry)
2109        (setf (bignum-ref res i) result-digit
2110              carry carry-out)))))
2111
2112(defun bignum-negate-to-pointer (big len res)
2113  (declare (fixnum len))
2114  (let* ((carry 1))
2115    (do* ((i 0 (1+ i))
2116          (j 0 (+ j 4)))
2117         ((= i len) carry)
2118      (declare (fixnum i))
2119      (multiple-value-bind (result-digit carry-out)
2120          (%add-with-carry (%lognot (bignum-ref big i)) 0 carry)
2121        (setf (%get-unsigned-long res j) result-digit
2122              carry carry-out)))))
2123 
2124
2125(defun %bignum-count-trailing-zero-bits (bignum)
2126  (let* ((count 0))
2127    (dotimes (i (%bignum-length bignum))
2128      (let* ((digit (bignum-ref bignum i)))
2129        (declare (type bignum-element-type digit))
2130        (if (zerop digit)
2131          (incf count 32)
2132          (progn
2133            (dotimes (bit 32)
2134              (declare (type (mod 32) bit))
2135              (if (logbitp bit digit)
2136                (return)
2137                (incf count)))
2138            (return)))))
2139    count))
2140                 
2141
2142(defun one-bignum-factor-of-two (a) 
2143  (declare (type bignum-type a))
2144  (let ((len (%bignum-length a)))
2145    (declare (fixnum len))
2146    (dotimes (i len)
2147      (let* ((x (bignum-ref a i)))
2148        (declare (fixnum x))
2149        (unless (zerop x)
2150          (return (+ (ash i 5)
2151                     (dotimes (j 32)
2152                       (if (logbitp j x)
2153                         (return j))))))))))
2154
2155
2156(defun %bignum-random (number state)
2157  (let* ((ndigits (%bignum-length number))
2158         (sign-index (1- ndigits)))
2159    (declare (fixnum ndigits sign-index))
2160    (with-bignum-buffers ((bignum ndigits))
2161      (dotimes (i sign-index)
2162        (setf (bignum-ref bignum i) (random (expt 2 digit-size) state)))
2163      (setf (bignum-ref bignum sign-index)
2164            (logand #x7fffffff (the (unsigned-byte 32)
2165                                 (random (expt 2 (1- digit-size)) state))))
2166      (let* ((result (mod bignum number)))
2167        (if (eq result bignum)
2168          (copy-uvector bignum)
2169          result)))))
2170
2171
2172
2173(defun logbitp (index integer)
2174  "Predicate returns T if bit index of integer is a 1."
2175  (number-case index
2176    (fixnum
2177     (if (minusp (the fixnum index))(report-bad-arg index '(integer 0))))
2178    (bignum
2179     ;; assuming bignum cant have more than most-positive-fixnum bits
2180     ;; (2 expt 24 longs)
2181     (if (bignum-minusp index)(report-bad-arg index '(integer 0)))
2182     ;; should error if integer isn't
2183     (return-from logbitp (minusp (require-type integer 'integer)))))
2184  (number-case integer
2185    (fixnum
2186     (if (%i< index (- target::nbits-in-word target::fixnumshift))
2187       (%ilogbitp index integer)
2188       (minusp (the fixnum integer))))
2189    (bignum
2190     (let ((bidx (%iasr 5 index))
2191           (bbit (%ilogand index 31)))
2192       (declare (fixnum bidx bbit))
2193       (if (>= bidx (%bignum-length integer))
2194         (bignum-minusp integer)
2195         (logbitp bbit (bignum-ref integer bidx)))))))
2196
2197) ; #+64-bit-target
Note: See TracBrowser for help on using the repository browser.