source: trunk/source/level-0/l0-bignum32.lisp @ 13066

Last change on this file since 13066 was 13066, checked in by rme, 10 years ago

Change "OpenMCL" to "Clozure CL" in comments and docstrings.

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