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

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

Update copyright notices.

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