source: branches/ia32/level-0/l0-bignum32.lisp @ 9370

Last change on this file since 9370 was 9370, checked in by rme, 13 years ago

most-positive-fixnum => target::target-most-positive-fixnum

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