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

Last change on this file since 10157 was 10157, checked in by rme, 12 years ago

On x8632, don't use fancy bignum multiplication algorithm. (The LAP
functions needed to use it aren't yet implemented.)

  • 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                        #+x8632-target
684                        nil)
685                 (let* ((ubytes (* len-a 4))
686                        (vbytes (* len-b 4))
687                        (rbytes (* len-res 4)))
688                   (declare (fixnum ubytes vbytes rbytes))
689                   (%stack-block ((uptr ubytes)
690                                  (vptr vbytes)
691                                  (rptr rbytes))
692                     (let* ((up (macptr->fixnum uptr))
693                            (vp (macptr->fixnum vptr))
694                            (rp (macptr->fixnum rptr)))
695                       (declare (fixnum up vp rp))
696                       (%copy-ivector-to-ptr a 0 uptr 0 ubytes)
697                       (if (eq a b)     ; maybe try eql ..
698                         (mpn-mul rp up len-a up len-a)
699                         (progn
700                           (%copy-ivector-to-ptr b 0 vptr 0 vbytes)
701                           (if (< len-a len-b)
702                             (mpn-mul rp vp len-b up len-a)
703                             (mpn-mul rp up len-a vp len-b)))))
704                     (%copy-ptr-to-ivector rptr 0 res 0 rbytes)))
705                 (dotimes (i len-a)
706                   (declare (type bignum-index i))
707                   (%multiply-and-add-harder-loop-2 a b res i len-b)))
708                 res)))
709      (let* ((res (with-negated-bignum-buffers a b multiply-unsigned-bignums)))
710        (if signs-differ (negate-bignum-in-place res))
711        (%normalize-bignum-macro res)))))
712
713
714(defun multiply-bignum-and-fixnum (bignum fixnum)
715  (declare (type bignum-type bignum) (fixnum fixnum))
716  (let* ((bignum-len (%bignum-length bignum))
717         (bignum-plus-p (bignum-plusp bignum))
718         (fixnum-plus-p (not (minusp fixnum)))
719         (negate-res (neq bignum-plus-p fixnum-plus-p)))
720    (declare (type bignum-type bignum)
721             (type bignum-index bignum-len))
722    (flet ((do-it (bignum fixnum  negate-res)
723             (let* ((bignum-len (%bignum-length bignum))
724                    (result (%allocate-bignum (the fixnum (1+ bignum-len)))))
725               (declare (type bignum-type bignum)
726                        (type bignum-index bignum-len))
727               (with-small-bignum-buffers ((carry-digit)
728                                           (result-digit))
729                 (dotimes (i bignum-len (%set-digit result bignum-len carry-digit))
730                   (%set-digit result i
731                               (%multiply-and-add result-digit carry-digit bignum i fixnum))))
732               (when negate-res
733                 (negate-bignum-in-place result))
734               (%normalize-bignum-macro result ))))
735      (declare (dynamic-extent #'do-it))
736      (if bignum-plus-p
737        (do-it bignum (if fixnum-plus-p fixnum (- fixnum))  negate-res)
738        (with-bignum-buffers ((b1 (the fixnum (1+ bignum-len))))
739          (negate-bignum bignum nil b1)
740          (do-it b1 (if fixnum-plus-p fixnum (- fixnum))  negate-res))))))
741
742;; assume we already know result won't fit in a fixnum
743;; only caller is fixnum-*-2
744;;
745
746(defun multiply-fixnums (a b)
747  (declare (fixnum a b))
748  (* a b))
749
750
751;;;; GCD.
752
753
754;;; Both args are > 0.
755(defun bignum-fixnum-gcd (bignum fixnum)
756  (let* ((rem (bignum-truncate-by-fixnum-no-quo bignum fixnum)))
757    (declare (fixnum rem))
758    (if (zerop rem)
759      fixnum
760      (%fixnum-gcd rem fixnum))))
761
762
763
764;;; NEGATE-BIGNUM -- Public.
765;;;
766;;; Fully-normalize is an internal optional.  It cause this to always return
767;;; a bignum, without any extraneous digits, and it never returns a fixnum.
768;;;
769(defun negate-bignum (x &optional (fully-normalize t) res)
770  (declare (type bignum-type x))
771  (let* ((len-x (%bignum-length x))
772         (len-res (1+ len-x))
773         (minusp (bignum-minusp x)))
774    (declare (type bignum-index len-x len-res))
775    (if (not res) (setq res (%allocate-bignum len-res))) ;Test len-res for range?
776    (let ((carry (bignum-negate-loop-really x len-x res)))  ; i think carry is always 0
777      (if (eq carry 0)
778        (if minusp (%bignum-set res len-x 0 0)(%bignum-set res len-x #xffff #xffff))
779        (digit-bind (h l)
780                    (if minusp 
781                      (%add-the-carry 0 0 carry)
782                      (%add-the-carry #xffff #xffff carry))
783                   
784          (%bignum-set res len-x h l))))
785    (if fully-normalize
786      (%normalize-bignum-macro res)
787      (%mostly-normalize-bignum-macro res))))
788
789;;; NEGATE-BIGNUM-IN-PLACE -- Internal.
790;;;
791;;; This assumes bignum is positive; that is, the result of negating it will
792;;; stay in the provided allocated bignum.
793;;;
794(defun negate-bignum-in-place (bignum)
795  (bignum-negate-loop-really bignum (%bignum-length bignum) bignum)
796  bignum)
797
798
799 
800
801(defun copy-bignum (bignum)
802  (let ((res (%allocate-bignum (%bignum-length bignum))))
803    (bignum-replace res bignum)
804    res))
805
806
807
808;;; BIGNUM-ASHIFT-RIGHT -- Public.
809;;;
810;;; First compute the number of whole digits to shift, shifting them by
811;;; skipping them when we start to pick up bits, and the number of bits to
812;;; shift the remaining digits into place.  If the number of digits is greater
813;;; than the length of the bignum, then the result is either 0 or -1.  If we
814;;; shift on a digit boundary (that is, n-bits is zero), then we just copy
815;;; digits.  The last branch handles the general case which uses a macro that a
816;;; couple other routines use.  The fifth argument to the macro references
817;;; locals established by the macro.
818;;;
819
820
821(defun bignum-ashift-right (bignum x)
822  (declare (type bignum-type bignum)
823           (fixnum x))
824  (let ((bignum-len (%bignum-length bignum)))
825    (declare (type bignum-index bignum-len))
826    (multiple-value-bind (digits n-bits) (truncate x digit-size)
827      (declare (type bignum-index digits)(fixnum n-bits))
828      (cond
829       ((>= digits bignum-len)
830        (if (bignum-plusp bignum) 0 -1))
831       ((eql 0 n-bits)
832        (bignum-ashift-right-digits bignum digits))
833       (t
834        (let* ((res-len (- bignum-len digits))
835               (res (%allocate-bignum res-len))
836               (len-1 (1- res-len)))
837          (declare (fixnum res-len len-1))
838          (bignum-shift-right-loop-1 n-bits res bignum len-1 digits)         
839          (%normalize-bignum-macro res )))))))
840
841                               
842
843
844
845;;; BIGNUM-ASHIFT-RIGHT-DIGITS -- Internal.
846;;;
847(defun bignum-ashift-right-digits (bignum digits)
848  (declare (type bignum-type bignum)
849           (type bignum-index digits))
850  (let* ((res-len (- (%bignum-length bignum) digits))
851         (res (%allocate-bignum res-len)))
852    (declare (type bignum-index res-len)
853             (type bignum-type res))
854    (bignum-replace res bignum :start2 digits)
855    (%normalize-bignum-macro res)))
856
857
858;;; BIGNUM-BUFFER-ASHIFT-RIGHT -- Internal.
859;;;
860;;; GCD uses this for an in-place shifting operation.  This is different enough
861;;; from BIGNUM-ASHIFT-RIGHT that it isn't worth folding the bodies into a
862;;; macro, but they share the basic algorithm.  This routine foregoes a first
863;;; test for digits being greater than or equal to bignum-len since that will
864;;; never happen for its uses in GCD.  We did fold the last branch into a macro
865;;; since it was duplicated a few times, and the fifth argument to it
866;;; references locals established by the macro.
867;;;
868#|
869(defun bignum-buffer-ashift-right (bignum bignum-len x)
870  (declare (type bignum-index bignum-len) (fixnum x))
871  (multiple-value-bind (digits n-bits)
872                       (truncate x digit-size)
873    (declare (type bignum-index digits))
874    (cond
875     ((zerop n-bits)
876      (let ((new-end (- bignum-len digits)))
877        (bignum-replace bignum bignum :end1 new-end :start2 digits
878                        :end2 bignum-len)
879        (%normalize-bignum-buffer bignum new-end)))
880     (t
881      (shift-right-unaligned bignum digits n-bits (- bignum-len digits)
882                             ((= j res-len-1)
883                              (digit-bind (h l) (%bignum-ref bignum i)
884                                (digit-set (h l) (%ashr h l n-bits))
885                                (%bignum-set bignum j h l))
886                              (%normalize-bignum-buffer bignum res-len)))))))
887|#
888#|
889(defun bignum-buffer-ashift-right (bignum bignum-len x)
890  (declare (type bignum-index bignum-len) (fixnum x))
891  (multiple-value-bind (digits n-bits) (truncate x digit-size)
892    (declare (type bignum-index digits)(fixnum n-bits))
893    (macrolet ((clear-high-digits ()
894                 `(do* ((i (1- (the fixnum (%bignum-length bignum))) (1- i))
895                        (j digits (1- j)))
896                       ((= 0 j))
897                    (declare (fixnum i j))
898                    (%bignum-set bignum i 0 0))))
899      (cond
900       ((zerop n-bits)
901        (let* ((new-end (- bignum-len digits)))
902          (declare (fixnum new-end))
903          (bignum-replace bignum bignum :end1 new-end :start2 digits
904                          :end2 bignum-len)
905          (clear-high-digits)
906          (%normalize-bignum-buffer bignum new-end)))
907       (t
908        (let* ((res-len (- bignum-len digits))
909               (len-1 (1- res-len)))
910          (declare (fixnum res-len len-1))
911          (bignum-shift-right-loop-1 n-bits bignum bignum len-1 digits)
912          ; clear the old high order digits - assume always positive
913          ; (when (neq 0 digits)(push digits poof))
914          (clear-high-digits)
915          (%normalize-bignum-buffer bignum res-len)))))))
916|#
917
918 
919
920;;; BIGNUM-ASHIFT-LEFT -- Public.
921;;;
922;;; This handles shifting a bignum buffer to provide fresh bignum data for some
923;;; internal routines.  We know bignum is safe when called with bignum-len.
924;;; First we compute the number of whole digits to shift, shifting them
925;;; starting to store farther along the result bignum.  If we shift on a digit
926;;; boundary (that is, n-bits is zero), then we just copy digits.  The last
927;;; branch handles the general case.
928;;;
929(defun bignum-ashift-left (bignum x &optional bignum-len)
930  (declare (type bignum-type bignum)
931           (fixnum x)
932           (type (or null bignum-index) bignum-len))
933  (multiple-value-bind (digits n-bits)
934                       (truncate x digit-size)
935    (declare (fixnum digits n-bits))
936    (let* ((bignum-len (or bignum-len (%bignum-length bignum)))
937           (res-len (+ digits bignum-len 1)))
938      (declare (fixnum bignum-len res-len))
939      (when (> res-len maximum-bignum-length)
940        (error "Can't represent result of left shift."))
941      (if (zerop n-bits)
942        (bignum-ashift-left-digits bignum bignum-len digits)
943        (bignum-ashift-left-unaligned bignum digits n-bits res-len)))))
944
945;;; BIGNUM-ASHIFT-LEFT-DIGITS -- Internal.
946;;;
947(defun bignum-ashift-left-digits (bignum bignum-len digits)
948  (declare (type bignum-index bignum-len digits))
949  (let* ((res-len (+ bignum-len digits))
950         (res (%allocate-bignum res-len)))
951    (declare (type bignum-index res-len))
952    (bignum-replace res bignum :start1 digits :end1 res-len :end2 bignum-len
953                    :from-end t)
954    res))
955
956;;; BIGNUM-ASHIFT-LEFT-UNALIGNED -- Internal.
957;;;
958;;; BIGNUM-TRUNCATE uses this to store into a bignum buffer by supplying res.
959;;; When res comes in non-nil, then this foregoes allocating a result, and it
960;;; normalizes the buffer instead of the would-be allocated result.
961;;;
962;;; We start storing into one digit higher than digits, storing a whole result
963;;; digit from parts of two contiguous digits from bignum.  When the loop
964;;; finishes, we store the remaining bits from bignum's first digit in the
965;;; first non-zero result digit, digits.  We also grab some left over high
966;;; bits from the last digit of bignum.
967;;;
968
969(defun bignum-ashift-left-unaligned (bignum digits n-bits res-len
970                                              &optional (res nil resp))
971  (declare (type bignum-index digits res-len)
972           (type (mod #.digit-size) n-bits))
973  (let* (;(remaining-bits (- digit-size n-bits))
974         (res-len-1 (1- res-len))
975         (res (or res (%allocate-bignum res-len))))
976    (declare (type bignum-index res-len res-len-1))
977    (bignum-shift-left-loop n-bits res bignum res-len-1 (the fixnum (1+ digits)))
978    ; if resp provided we don't care about returned value
979    (if (not resp) (%normalize-bignum-macro res))))
980
981
982
983
984;;;; Relational operators.
985
986
987
988;;; BIGNUM-COMPARE -- Public.
989;;;
990;;; This compares two bignums returning -1, 0, or 1, depending on whether a
991;;; is less than, equal to, or greater than b.
992;;;
993;(proclaim '(function bignum-compare (bignum bignum) (integer -1 1)))
994(defun bignum-compare (a b)
995  (declare (type bignum-type a b))
996  (let* ((a-plusp (bignum-plusp a))
997         (b-plusp (bignum-plusp b)))
998    (if (eq a-plusp b-plusp)
999      (let* ((len-a (%bignum-length a))
1000             (len-b (%bignum-length b)))
1001        (declare (type bignum-index len-a len-b))
1002        (cond ((= len-a len-b)
1003               (do* ((i (1- len-a) (1- i)))
1004                    ((zerop i) (%compare-digits a b 0))
1005                 (declare (fixnum i))
1006                 (let* ((signum (%compare-digits a b i)))
1007                   (declare (fixnum signum))
1008                   (unless (zerop signum)
1009                     (return signum)))))
1010              ((> len-a len-b)
1011               (if a-plusp 1 -1))
1012              (t (if a-plusp -1 1))))
1013      (if a-plusp 1 -1))))
1014
1015
1016
1017
1018
1019
1020;;;; Integer length and logcount
1021
1022
1023(defun bignum-integer-length (big)
1024  (the fixnum (- (the fixnum (ash (the fixnum (%bignum-length big)) 5))
1025                 (the fixnum (%bignum-sign-bits big)))))
1026
1027; (not (zerop (logand integer1 integer2)
1028
1029(defun bignum-logtest (num1 num2)
1030  (let* ((length1 (%bignum-length num1))
1031         (length2 (%bignum-length num2))
1032         (n1-minusp (bignum-minusp num1))
1033         (n2-minusp (bignum-minusp num2)))
1034    (declare (fixnum length1 length2))
1035    (if (and n1-minusp n2-minusp) ; both neg, get out quick
1036      T       
1037      (let ((val (bignum-logtest-loop (min length1 length2) num1 num2)))
1038                 #|(do* ((index 0 (1+ index)))
1039                      ((= index (min length1 length2)) nil)
1040                   ; maybe better to start from high end of shorter?
1041                   (multiple-value-bind (hi1 lo1)(%bignum-ref num1 index)
1042                     (multiple-value-bind (hi2 lo2)(%bignum-ref num2 index)
1043                       (when (or (not (zerop (%ilogand hi1 hi2)))
1044                                 (not (zerop (%ilogand lo1 lo2))))
1045                         (return t)))))))|#
1046        (or val
1047            (when (not (eql length1 length2)) ; lengths same => value nil
1048              (if (< length1 length2)
1049                n1-minusp
1050                n2-minusp)))))))
1051
1052
1053
1054(defun logtest-fix-big (fix big)
1055  (declare (fixnum fix))
1056  (if (eql 0 (the fixnum fix))
1057    nil
1058    (if (> (the fixnum fix) 0) 
1059      (let ()
1060        (multiple-value-bind (hi lo)(%bignum-ref big 0)
1061          (declare (fixnum hi lo))
1062          (or (not (zerop (logand fix lo)))
1063              (not (zerop (logand (ash fix (- 16)) hi))))))
1064      t)))
1065
1066
1067(defun bignum-logcount (bignum)
1068  (declare (type bignum-type bignum))
1069  (let* ((length (%bignum-length bignum))
1070         (plusp (bignum-plusp bignum))
1071         (result 0))
1072    (declare (type bignum-index length)
1073             (fixnum result))
1074    (if plusp
1075      (dotimes (index length result)
1076        (incf result (the fixnum (%logcount bignum index))))
1077      (dotimes (index length result)
1078        (incf result (the fixnum (%logcount-complement bignum index)))))))
1079
1080
1081;;;; Logical operations.
1082
1083;;; NOT.
1084;;;
1085
1086;;; BIGNUM-LOGICAL-NOT -- Public.
1087;;;
1088(defun bignum-logical-not (a)
1089  (declare (type bignum-type a))
1090  (let* ((len (%bignum-length a))
1091         (res (%allocate-bignum len)))
1092    (declare (type bignum-index len))
1093    (dotimes (i len res)
1094      (%bignum-lognot i a res))))
1095
1096
1097
1098
1099;;; AND.
1100;;;
1101
1102;;; BIGNUM-LOGICAL-AND -- Public.
1103;;;
1104(defun bignum-logical-and (a b)
1105  (declare (type bignum-type a b))
1106  (let* ((len-a (%bignum-length a))
1107         (len-b (%bignum-length b))
1108         (a-plusp (bignum-plusp a))
1109         (b-plusp (bignum-plusp b)))
1110    (declare (type bignum-index len-a len-b))
1111    (cond
1112      ((< len-a len-b)
1113       (if a-plusp
1114         (logand-shorter-positive a len-a b (%allocate-bignum len-a))
1115         (logand-shorter-negative a len-a b len-b (%allocate-bignum len-b))))
1116      ((< len-b len-a)
1117       (if b-plusp
1118         (logand-shorter-positive b len-b a (%allocate-bignum len-b))
1119         (logand-shorter-negative b len-b a len-a (%allocate-bignum len-a))))
1120      (t (logand-shorter-positive a len-a b (%allocate-bignum len-a))))))
1121
1122;;; LOGAND-SHORTER-POSITIVE -- Internal.
1123;;;
1124;;; This takes a shorter bignum, a and len-a, that is positive.  Because this
1125;;; is AND, we don't care about any bits longer than a's since its infinite 0
1126;;; sign bits will mask the other bits out of b.  The result is len-a big.
1127;;;
1128(defun logand-shorter-positive (a len-a b res)
1129  (declare (type bignum-type a b res)
1130           (type bignum-index len-a))
1131  (dotimes (i len-a)
1132    (%bignum-logand i a b res))
1133  (%normalize-bignum-macro res))
1134
1135;;; LOGAND-SHORTER-NEGATIVE -- Internal.
1136;;;
1137;;; This takes a shorter bignum, a and len-a, that is negative.  Because this
1138;;; is AND, we just copy any bits longer than a's since its infinite 1 sign
1139;;; bits will include any bits from b.  The result is len-b big.
1140;;;
1141(defun logand-shorter-negative (a len-a b len-b res)
1142  (declare (type bignum-type a b res)
1143           (type bignum-index len-a len-b))
1144  (dotimes (i len-a)
1145    (%bignum-logand i a b res))
1146  (bignum-replace res b :start1 len-a :start2 len-a :end1 len-b :end2 len-b)
1147  (%normalize-bignum-macro res))
1148
1149
1150
1151;;;
1152;;;
1153;;; bignum-logandc2
1154
1155(defun bignum-logandc2 (a b)
1156  (declare (type bignum-type a b))
1157  (let* ((len-a (%bignum-length a))
1158         (len-b (%bignum-length b))
1159         (a-plusp (bignum-plusp a))
1160         (b-plusp (bignum-plusp b)))
1161    (declare (type bignum-index len-a len-b))
1162    (cond
1163     ((< len-a len-b)
1164      (logandc2-shorter-any a len-a b len-b (if a-plusp (%allocate-bignum len-a) (%allocate-bignum len-b))))
1165     ((< len-b len-a) ; b shorter
1166      (logandc1-shorter-any b len-b a len-a (if b-plusp (%allocate-bignum len-a)(%allocate-bignum len-b))))
1167     (t (logandc2-shorter-any a len-a b len-b (%allocate-bignum len-a))))))
1168
1169(defun logandc2-shorter-any (a len-a b len-b res)
1170  (declare (type bignum-type a b res)
1171           (type bignum-index len-a len-b))
1172  (dotimes (i len-a)
1173    (%bignum-logandc2 i a b res))
1174  (if (bignum-minusp a)
1175    (do ((i len-a (1+ i)))
1176          ((= i len-b))
1177        (declare (type bignum-index i))
1178        (digit-bind (h l) (%bignum-ref b i)
1179          (%bignum-set res i (%ilognot h) (%ilognot l)))))
1180  (%normalize-bignum-macro res))
1181
1182
1183
1184(defun logandc1-shorter-any (a len-a b len-b res)
1185  (declare (type bignum-type a b res)
1186           (type bignum-index len-a len-b))
1187  (dotimes (i len-a)
1188    (%bignum-logandc1 i a b res))
1189  (if (bignum-plusp a)
1190    (if (neq len-a len-b)
1191      (bignum-replace res b :start1 len-a :start2 len-a :end1 len-b :end2 len-b)))
1192  (%normalize-bignum-macro res))
1193
1194
1195
1196(defun fix-big-logand (fix big)
1197  (let* ((len-b (%bignum-length big))
1198         (res (if (< fix 0)(%allocate-bignum len-b))))
1199    (declare (fixnum fix len-b))       
1200    (let ((val (fix-digit-logand fix big res)))
1201      (if res
1202        (progn
1203          (bignum-replace res big :start1 1 :start2 1 :end1 len-b :end2 len-b)
1204          (%normalize-bignum-macro res))
1205        val))))
1206 
1207
1208(defun fix-big-logandc2 (fix big)
1209  (let* ((len-b (%bignum-length big))
1210         (res (if (< fix 0)(%allocate-bignum len-b))))
1211    (declare (fixnum fix len-b))       
1212    (let ((val (fix-digit-logandc2 fix big res)))
1213      (if res
1214        (progn
1215          (do ((i 1 (1+ i)))
1216              ((= i len-b))
1217            (declare (type bignum-index i))
1218            (digit-lognot-move i big res))
1219          (%normalize-bignum-macro res))
1220        val))))
1221
1222(defun fix-big-logandc1 (fix big)
1223  (let* ((len-b (%bignum-length big))
1224         (res (if (>= fix 0)(%allocate-bignum len-b))))
1225    (declare (fixnum fix len-b))       
1226    (let ((val (fix-digit-logandc1 fix big res)))
1227      (if res
1228        (progn 
1229          (bignum-replace res big :start1 1 :start2 1 :end1 len-b :end2 len-b)
1230          (%normalize-bignum-macro res))
1231        val))))
1232
1233
1234
1235
1236
1237
1238
1239;;; IOR.
1240;;;
1241
1242;;; BIGNUM-LOGICAL-IOR -- Public.
1243;;;
1244(defun bignum-logical-ior (a b)
1245  (declare (type bignum-type a b))
1246  (let* ((len-a (%bignum-length a))
1247         (len-b (%bignum-length b))
1248         (a-plusp (bignum-plusp a))
1249         (b-plusp (bignum-plusp b)))
1250    (declare (type bignum-index len-a len-b))
1251    (cond
1252     ((< len-a len-b)
1253      (if a-plusp
1254          (logior-shorter-positive a len-a b len-b (%allocate-bignum len-b))
1255          (logior-shorter-negative a len-a b len-b (%allocate-bignum len-b))))
1256     ((< len-b len-a)
1257      (if b-plusp
1258          (logior-shorter-positive b len-b a len-a (%allocate-bignum len-a))
1259          (logior-shorter-negative b len-b a len-a (%allocate-bignum len-a))))
1260     (t (logior-shorter-positive a len-a b len-b (%allocate-bignum len-a))))))
1261
1262;;; LOGIOR-SHORTER-POSITIVE -- Internal.
1263;;;
1264;;; This takes a shorter bignum, a and len-a, that is positive.  Because this
1265;;; is IOR, we don't care about any bits longer than a's since its infinite
1266;;; 0 sign bits will mask the other bits out of b out to len-b.  The result
1267;;; is len-b long.
1268;;;
1269(defun logior-shorter-positive (a len-a b len-b res)
1270  (declare (type bignum-type a b res)
1271           (type bignum-index len-a len-b))
1272  (dotimes (i len-a)
1273    (%bignum-logior i a b res))
1274  (if (not (eql len-a len-b))
1275    (bignum-replace res b :start1 len-a :start2 len-a :end1 len-b :end2 len-b))
1276  (%normalize-bignum-macro res))
1277
1278;;; LOGIOR-SHORTER-NEGATIVE -- Internal.
1279;;;
1280;;; This takes a shorter bignum, a and len-a, that is negative.  Because this
1281;;; is IOR, we just copy any bits longer than a's since its infinite 1 sign
1282;;; bits will include any bits from b.  The result is len-b long.
1283;;;
1284(defun logior-shorter-negative (a len-a b len-b res)
1285  (declare (type bignum-type a b res)
1286           (type bignum-index len-a len-b))
1287  (dotimes (i len-a)
1288    (%bignum-logior i a b res))
1289  ; silly to propagate sign and then normalize it away
1290  ; but may need to do at least once - but we are only normalizing from len-a?
1291  ; ah but the sign needs to be correct
1292  (do ((i len-a (1+ i)))
1293      ((= i len-b))
1294    (declare (type bignum-index i))
1295    (%bignum-set res i #xffff #xffff))
1296  (%normalize-bignum-macro res))
1297
1298
1299
1300
1301;;; XOR.
1302;;;
1303
1304;;; BIGNUM-LOGICAL-XOR -- Public.
1305;;;
1306(defun bignum-logical-xor (a b)
1307  (declare (type bignum-type a b))
1308  (let ((len-a (%bignum-length a))
1309        (len-b (%bignum-length b)))
1310    (declare (type bignum-index len-a len-b))
1311    (if (< len-a len-b)
1312        (bignum-logical-xor-aux a len-a b len-b (%allocate-bignum len-b))
1313        (bignum-logical-xor-aux b len-b a len-a (%allocate-bignum len-a)))))
1314
1315;;; BIGNUM-LOGICAL-XOR-AUX -- Internal.
1316;;;
1317;;; This takes the the shorter of two bignums in a and len-a.  Res is len-b
1318;;; long.  Do the XOR.
1319;;;
1320(defun bignum-logical-xor-aux (a len-a b len-b res)
1321  (declare (type bignum-type a b res)
1322           (type bignum-index len-a len-b))
1323  (dotimes (i len-a)
1324    (%bignum-logxor i a b res))
1325  (unless (= len-a len-b)
1326    (let ((sign (if (bignum-minusp a) #xffff 0)))
1327      (do ((i len-a (1+ i)))
1328          ((= i len-b))
1329        (declare (type bignum-index i))
1330        (digit-bind (h l) (%bignum-ref b i)
1331          (%bignum-set res i (%ilogxor sign h)(%ilogxor sign l))))))
1332  (%normalize-bignum-macro res))
1333
1334
1335
1336
1337
1338;;;; LDB (load byte)
1339
1340; [slh] 'twas all commented out - thank gawd
1341
1342
1343;;;; TRUNCATE.
1344
1345;;; This is the original sketch of the algorithm from which I implemented this
1346;;; TRUNCATE, assuming both operands are bignums.  I should modify this to work
1347;;; with the documentation on my functions, as a general introduction.  I've
1348;;; left this here just in case someone needs it in the future.  Don't look
1349;;; at this unless reading the functions' comments leaves you at a loss.
1350;;; Remember this comes from Knuth, so the book might give you the right general
1351;;; overview.
1352;;;
1353;;;
1354;;; (truncate x y):
1355;;;
1356;;; If X's magnitude is less than Y's, then result is 0 with remainder X.
1357;;;
1358;;; Make x and y positive, copying x if it is already positive.
1359;;;
1360;;; Shift y left until there's a 1 in the 30'th bit (most significant, non-sign
1361;;;       digit)
1362;;;    Just do most sig digit to determine how much to shift whole number.
1363;;; Shift x this much too.
1364;;; Remember this initial shift count.
1365;;;
1366;;; Allocate q to be len-x minus len-y quantity plus 1.
1367;;;
1368;;; i = last digit of x.
1369;;; k = last digit of q.
1370;;;
1371;;; LOOP
1372;;;
1373;;; j = last digit of y.
1374;;;
1375;;; compute guess.
1376;;; if x[i] = y[j] then g = #xFFFFFFFF
1377;;; else g = x[i]x[i-1]/y[j].
1378;;;
1379;;; check guess.
1380;;; %UNSIGNED-MULTIPLY returns b and c defined below.
1381;;;    a = x[i-1] - (logand (* g y[j]) #xFFFFFFFF).
1382;;;       Use %UNSIGNED-MULTIPLY taking low-order result.
1383;;;    b = (logand (ash (* g y[j-1]) -32) #xFFFFFFFF).
1384;;;    c = (logand (* g y[j-1]) #xFFFFFFFF).
1385;;; if a < b, okay.
1386;;; if a > b, guess is too high
1387;;;    g = g - 1; go back to "check guess".
1388;;; if a = b and c > x[i-2], guess is too high
1389;;;    g = g - 1; go back to "check guess".
1390;;; GUESS IS 32-BIT NUMBER, SO USE THING TO KEEP IN SPECIAL REGISTER
1391;;; SAME FOR A, B, AND C.
1392;;;
1393;;; Subtract g * y from x[i - len-y+1]..x[i].  See paper for doing this in step.
1394;;; If x[i] < 0, guess is fucked.
1395;;;    negative g, then add 1
1396;;;    zero or positive g, then subtract 1
1397;;; AND add y back into x[len-y+1..i].
1398;;;
1399;;; q[k] = g.
1400;;; i = i - 1.
1401;;; k = k - 1.
1402;;;
1403;;; If k>=0, goto LOOP.
1404;;;
1405;;;
1406;;; Now quotient is good, but remainder is not.
1407;;; Shift x right by saved initial left shifting count.
1408;;;
1409;;; Check quotient and remainder signs.
1410;;; x pos y pos --> q pos r pos
1411;;; x pos y neg --> q neg r pos
1412;;; x neg y pos --> q neg r neg
1413;;; x neg y neg --> q pos r neg
1414;;;
1415;;; Normalize quotient and remainder.  Cons result if necessary.
1416;;;
1417
1418
1419
1420;;; These are used by BIGNUM-TRUNCATE and friends in the general case.
1421;;;
1422(defvar *truncate-x* nil)
1423(defvar *truncate-y* nil)
1424
1425;;; BIGNUM-TRUNCATE -- Public.
1426;;;
1427;;; This divides x by y returning the quotient and remainder.  In the general
1428;;; case, we shift y to setup for the algorithm, and we use two buffers to save
1429;;; consing intermediate values.  X gets destructively modified to become the
1430;;; remainder, and we have to shift it to account for the initial Y shift.
1431;;; After we multiple bind q and r, we first fix up the signs and then return
1432;;; the normalized results.
1433;;;
1434
1435
1436(defun bignum-truncate (x1 y1 &optional no-rem)
1437  (declare (type bignum-type x1 y1))
1438  (let* ((x-plusp (bignum-plusp x1))
1439         (y-plusp (bignum-plusp y1)))
1440    (flet 
1441      ((do-it (x y) 
1442         (let* ((len-x (%bignum-length x))
1443                (len-y (%bignum-length y)))
1444           (declare (fixnum len-x len-y))
1445           
1446           (let ((c (bignum-compare y x)))
1447             (cond 
1448              ((eql c 1)  ; >
1449               (return-from bignum-truncate (values 0 x1)))
1450              ((eql c 0)(values 1 0))  ; =  might as well since did compare anyway
1451              ((< len-y 2)
1452               (multiple-value-bind (q r)
1453                                    (bignum-truncate-single-digit x len-x y no-rem)
1454                 (values q
1455                         (unless no-rem
1456                           (cond (x-plusp r)
1457                                 ((typep r 'fixnum) (the fixnum (- (the fixnum r))))
1458                                 (t (negate-bignum-in-place r)
1459                                    (%normalize-bignum-macro r )))))))
1460              (t
1461               (let* ((len-x+1 (1+ len-x)))
1462                 (declare (fixnum len-x+1))
1463                 (with-bignum-buffers ((*truncate-x* len-x+1)
1464                                       (*truncate-y* (the fixnum (1+ len-y))))
1465                   (let ((y-shift (shift-y-for-truncate y)))
1466                     (shift-and-store-truncate-buffers x len-x y len-y y-shift)
1467                     (values (do-truncate len-x+1 len-y)
1468                             ;; DO-TRUNCATE must execute first.
1469                             (when (not no-rem)                               
1470                               (when (not (eql 0 y-shift))                                 
1471                                 (let* ((res-len-1 (1- len-y)))
1472                                   (declare (fixnum res-len-1))
1473                                   (bignum-shift-right-loop-1 y-shift *truncate-x* *truncate-x* res-len-1 0)))                               
1474                               (let ((the-res (%normalize-bignum-macro *truncate-x* )))
1475                                 (if (not (fixnump the-res))
1476                                   (if x-plusp (copy-bignum the-res) (negate-bignum the-res))
1477                                   (if x-plusp the-res (the fixnum (- (the fixnum the-res)))))
1478                                     ))))))))))))
1479      (multiple-value-bind (q r)(with-negated-bignum-buffers x1 y1 do-it)
1480        (let ((quotient (cond ((eq x-plusp y-plusp) q)
1481                              ((typep q 'fixnum) (the fixnum (- (the fixnum q))))
1482                              (t (negate-bignum-in-place q)
1483                                 (%normalize-bignum-macro q )))))
1484          (if no-rem
1485            quotient           
1486            (values quotient r)))))))
1487
1488(defun bignum-rem (x1 y1)
1489  (declare (type bignum-type x1 y1)) 
1490  (let* ((x-plusp (bignum-plusp x1)))
1491    (flet 
1492      ((do-it (x y) 
1493         (let* ((len-x (%bignum-length x))
1494                (len-y (%bignum-length y)))
1495           (declare (fixnum len-x len-y))           
1496           (let ((c (bignum-compare y x)))
1497             (cond 
1498              ((eql c 1) (return-from bignum-rem x1))
1499              ((eql c 0) 0)  ; =  might as well since did compare anyway
1500              ((< len-y 2)
1501               (let ((r (bignum-truncate-single-digit-no-quo x len-x y)))  ; phooey
1502                 (cond (x-plusp r)
1503                       ((typep r 'fixnum) (the fixnum (- (the fixnum r))))
1504                       (t (negate-bignum-in-place r)
1505                          (%normalize-bignum-macro r )))))
1506              (t
1507               (let* ((len-x+1 (1+ len-x)))
1508                 (declare (fixnum len-x+1))
1509                 (with-bignum-buffers ((*truncate-x* len-x+1)
1510                                       (*truncate-y* (the fixnum (1+ len-y))))
1511                   (let ((y-shift (shift-y-for-truncate y)))
1512                     (shift-and-store-truncate-buffers x len-x y len-y y-shift)
1513                     (do-truncate-no-quo len-x+1 len-y)
1514                     (when (not (eql 0 y-shift))                                 
1515                       (let* ((res-len-1 (1- len-y)))
1516                         (declare (fixnum res-len-1))
1517                         (bignum-shift-right-loop-1 y-shift *truncate-x* *truncate-x* res-len-1 0)))
1518                     (let ((the-res (%normalize-bignum-macro *truncate-x*)))
1519                       (if (not (fixnump the-res))
1520                         (if x-plusp (copy-bignum the-res) (negate-bignum the-res))
1521                         (if x-plusp the-res (the fixnum (- (the fixnum the-res)))))))))))))))
1522      (declare (dynamic-extent #'do-it))
1523      (with-negated-bignum-buffers x1 y1 do-it))))
1524
1525
1526
1527;;; BIGNUM-TRUNCATE-SINGLE-DIGIT -- Internal.
1528;;;
1529;;; This divides x by y when y is a single bignum digit.  BIGNUM-TRUNCATE fixes
1530;;; up the quotient and remainder with respect to sign and normalization.
1531;;;
1532;;; We don't have to worry about shifting y to make its most significant digit
1533;;; sufficiently large for %FLOOR to return 32-bit quantities for the q-digit
1534;;; and r-digit.  If y is a single digit bignum, it is already large enough
1535;;; for %FLOOR.  That is, it has some bits on pretty high in the digit.
1536;;;
1537;;; x is positive
1538(defun bignum-truncate-single-digit (x len-x y &optional no-rem)
1539  (declare (type bignum-index len-x))
1540  (let* ((maybe-q (%allocate-bignum 2))
1541         (q (if (<= len-x 2) maybe-q (%allocate-bignum len-x)))
1542         (r-h 0)
1543         (r-l 0))
1544    (declare (dynamic-extent maybe-q))
1545    (digit-bind (y-h y-l) (%bignum-ref y 0)
1546      (multiple-value-setq (r-h r-l)(%floor-loop-quo x q y-h y-l))     
1547      (if (eq q maybe-q)
1548        (progn 
1549          (setq q (%normalize-bignum-macro q))
1550          (if (not (fixnump q)) (setq q (copy-bignum q))))
1551        (setq q (%normalize-bignum-macro q )))
1552      ;; might as well make a fixnum if possible
1553      (if no-rem
1554        q
1555        (if (> (%digits-sign-bits r-h r-l)  target::fixnumshift)
1556          (values q (%ilogior (%ilsl 16 r-h) r-l))
1557          (let ((rem (%allocate-bignum 1)))
1558            (%bignum-set rem 0 r-h r-l)
1559            (values q rem)))))))
1560
1561;;; aka rem
1562(defun bignum-truncate-single-digit-no-quo (x len-x y)
1563  (declare (type bignum-index len-x))
1564  (declare (ignore len-x))
1565  (let (;(q (%allocate-bignum len-x))
1566        (r-h 0)
1567        (r-l 0))
1568    (progn
1569      (digit-bind (y-h y-l) (%bignum-ref y 0)
1570        (multiple-value-setq (r-h r-l)(%floor-loop-no-quo x y-h y-l))
1571        ; might as well make a fixnum if possible
1572        (if (> (%digits-sign-bits r-h r-l)  target::fixnumshift)
1573          (%ilogior (%ilsl 16 r-h) r-l)
1574          (let ((rem (%allocate-bignum 1)))
1575            (%bignum-set rem 0 r-h r-l)
1576            rem))))))
1577
1578;; so big deal - we save a one digit bignum for y
1579;; and bigger deal if x is negative - we copy or negate x, computing result destructively
1580;;  - thus avoiding making a negated x in addition to result
1581;;
1582(defun bignum-truncate-by-fixnum (x y)
1583  (declare (fixnum y))
1584  (when (eql y 0)(error (make-condition 'division-by-zero :operation 'truncate :operands (list x y))))
1585  (let* ((len-x (%bignum-length x))
1586         (x-minus (bignum-minusp x))
1587         (maybe-q (%allocate-bignum 3))
1588         (q (if x-minus
1589              (if (<= len-x 2)
1590                (dotimes (i 3 (negate-bignum-in-place maybe-q))
1591                  (if (< i len-x)
1592                    (multiple-value-bind (hi lo) (%bignum-ref x i)
1593                      (%bignum-set maybe-q i hi lo))
1594                    (%bignum-set maybe-q i 65535 65535)))
1595                (negate-bignum x))
1596              (if (<= len-x 2) ; this was broken if negative because bignum-replace just copies min len-a len-b digits
1597                (progn
1598                  (bignum-replace maybe-q x)               
1599                  maybe-q)
1600                (%allocate-bignum len-x))))      ;  q is new big or -x
1601         ;(len-q (%bignum-length q))
1602         (y-minus (minusp y))         
1603         (y (if y-minus (- y) y)))
1604    (declare (fixnum y))
1605    (declare (type bignum-index len-x))
1606    (declare (dynamic-extent maybe-q))
1607    (let* ((r-h 0)
1608           (r-l 0)
1609           (y-h (%ilogand #xffff (%iasr 16 y)))
1610           (y-l (%ilogand #xffff y)))
1611      (multiple-value-setq (r-h r-l)(%floor-loop-quo (if x-minus q x) q y-h y-l))     
1612      (let* ((r (%ilogior (%ilsl 16 r-h) r-l)))
1613        (declare (fixnum r))
1614        (when (neq x-minus y-minus)(negate-bignum-in-place q))
1615        (setq q (%normalize-bignum-macro q ))
1616        (values (if (eq q maybe-q) (copy-bignum q) q)
1617                (if x-minus (the fixnum (- r)) r))))))
1618
1619(defun bignum-truncate-by-fixnum-no-quo (x y)
1620  (declare (fixnum y))
1621  (when (eql y 0)(error (make-condition 'division-by-zero :operation 'truncate :operands (list x Y))))
1622  (let* ((len-x (%bignum-length x))
1623         (x-minus (bignum-minusp x))
1624         (y-minus (minusp y))         
1625         (y (if y-minus (- y) y)))
1626    (declare (fixnum y))
1627    (declare (type bignum-index len-x))
1628      (let* (;(LEN-Q (%BIGNUM-LENGTH Q))
1629             (r-h 0)
1630             (r-l 0)
1631             (y-h (%ilogand #xffff (%iasr 16 y)))
1632             (y-l (%ilogand #xffff y)))
1633        (if x-minus
1634          (with-bignum-buffers ((q (the fixnum (1+ len-x))))
1635            (negate-bignum x nil q)
1636            (multiple-value-setq (r-h r-l)(%floor-loop-no-quo q y-h y-l)))
1637          (multiple-value-setq (r-h r-l)(%floor-loop-no-quo x y-h y-l)))       
1638        (let* ((r (%ilogior (%ilsl 16 r-h) r-l)))
1639          (declare (fixnum r))
1640          (if x-minus (the fixnum (- r)) r)))))
1641
1642
1643;;; DO-TRUNCATE -- Internal.
1644;;;
1645;;; This divides *truncate-x* by *truncate-y*, and len-x and len-y tell us how
1646;;; much of the buffers we care about.  TRY-BIGNUM-TRUNCATE-GUESS modifies
1647;;; *truncate-x* on each interation, and this buffer becomes our remainder.
1648;;;
1649;;; *truncate-x* definitely has at least three digits, and it has one more than
1650;;; *truncate-y*.  This keeps i, i-1, i-2, and low-x-digit happy.  Thanks to
1651;;; SHIFT-AND-STORE-TRUNCATE-BUFFERS.
1652;;;
1653
1654
1655(defun do-truncate (len-x len-y)
1656  (declare (type bignum-index len-x len-y))
1657  (let* ((len-q (- len-x len-y))
1658         ;; Add one for extra sign digit in case high bit is on.
1659         (len-res (1+ len-q))
1660         (maybe-q (%allocate-bignum 2))         
1661         (q (if (<= len-res 2) maybe-q (%allocate-bignum len-res)))
1662         (k (1- len-q))
1663         (i (1- len-x))
1664         (low-x-digit (- i len-y)))
1665    (declare (type bignum-index len-q len-res k i  low-x-digit))
1666    (declare (dynamic-extent maybe-q))
1667    (loop
1668      (digit-bind (h l)
1669                  (digit-bind (guess-h guess-l)
1670                              (bignum-truncate-guess-2 *truncate-x* i *truncate-y* (the fixnum (1- len-y)))                                 
1671                    (try-bignum-truncate-guess guess-h guess-l len-y low-x-digit))
1672        (%bignum-set q k h l))
1673      (cond ((zerop k) (return))
1674            (t (decf k)
1675               (decf low-x-digit)
1676               (setq i (1- i)))))
1677    (if (eq q maybe-q)
1678      (progn 
1679        (setq q (%normalize-bignum-macro q))
1680        (if (fixnump q) q (copy-bignum q)))
1681      (%normalize-bignum-macro q))))
1682
1683(defun do-truncate-no-quo (len-x len-y)
1684  (declare (type bignum-index len-x len-y))
1685  (let* ((len-q (- len-x len-y))
1686         (k (1- len-q))
1687         (i (1- len-x))
1688         (low-x-digit (- i len-y)))
1689    (declare (type bignum-index len-q k i  low-x-digit))
1690    (loop
1691      (digit-bind (guess-h guess-l) (bignum-truncate-guess-2 *truncate-x* i *truncate-y* (the fixnum (1- len-y)))                                 
1692        (try-bignum-truncate-guess guess-h guess-l len-y low-x-digit)
1693        (cond ((zerop k) (return))
1694              (t (decf k)
1695                 (decf low-x-digit)
1696                 (setq i (1- i))))))
1697    nil))
1698
1699;;; TRY-BIGNUM-TRUNCATE-GUESS -- Internal.
1700;;;
1701;;; This takes a digit guess, multiplies it by *truncate-y* for a result one
1702;;; greater in length than len-y, and subtracts this result from *truncate-x*.
1703;;; Low-x-digit is the first digit of x to start the subtraction, and we know x
1704;;; is long enough to subtract a len-y plus one length bignum from it.  Next we
1705;;; check the result of the subtraction, and if the high digit in x became
1706;;; negative, then our guess was one too big.  In this case, return one less
1707;;; than guess passed in, and add one value of y back into x to account for
1708;;; subtracting one too many.  Knuth shows that the guess is wrong on the order
1709;;; of 3/b, where b is the base (2 to the digit-size power) -- pretty rarely.
1710;;;
1711
1712(defun try-bignum-truncate-guess (guess-h guess-l len-y low-x-digit)
1713  (declare (type bignum-index low-x-digit len-y))
1714
1715  (let ((carry-digit-h 0)
1716        (carry-digit-l 0)
1717        (borrow 1)
1718        (i low-x-digit))
1719    (declare (type bignum-index i)
1720             (fixnum borrow carry-digit-h carry-digit-l))
1721    ;; Multiply guess and divisor, subtracting from dividend simultaneously.
1722    (dotimes (j len-y)
1723      (multiple-value-bind (y-h y-l) (%bignum-ref *truncate-y* j)
1724        (multiple-value-bind (high-h high-l low-h low-l)
1725            (%multiply-and-add-1 guess-h
1726                               guess-l
1727                               y-h
1728                               y-l
1729                               carry-digit-h
1730                               carry-digit-l)
1731          (setq carry-digit-h high-h
1732                carry-digit-l high-l)
1733          (multiple-value-bind (tx-h tx-l) (%bignum-ref *truncate-x* i)
1734            (multiple-value-bind (x-h x-l temp-borrow)
1735                (%subtract-with-borrow-1 tx-h tx-l low-h low-l borrow)
1736              (%bignum-set *truncate-x* i x-h x-l)
1737              (setq borrow temp-borrow)))))
1738      (incf i))
1739    (multiple-value-bind (tx-h tx-l) (%bignum-ref *truncate-x* i)
1740      (multiple-value-bind (x-h x-l)
1741          (%subtract-with-borrow-1 tx-h tx-l carry-digit-h carry-digit-l borrow)
1742        (%bignum-set *truncate-x* i x-h x-l)))
1743    ;; See if guess is off by one, adding one Y back in if necessary.
1744
1745
1746    (cond ((%digit-0-or-plusp *truncate-x* i)
1747           (values guess-h guess-l))
1748          (t
1749           ;; If subtraction has negative result, add one divisor value back
1750           ;; in.  The guess was one too large in magnitude.
1751           ;; hmm - happens about 1.6% of the time
1752           (bignum-add-loop-+ low-x-digit *truncate-x* *truncate-y* len-y)
1753           (%subtract-one guess-h guess-l)
1754           ;(%subtract-with-borrow guess-h guess-l 0 1 1)
1755           ))))
1756
1757
1758
1759;;; BIGNUM-TRUNCATE-GUESS -- Internal.
1760;;;
1761;;; This returns a guess for the next division step.  Y1 is the highest y
1762;;; digit, and y2 is the second to highest y digit.  The x... variables are
1763;;; the three highest x digits for the next division step.
1764;;;
1765;;; From Knuth, our guess is either all ones or x-i and x-i-1 divided by y1,
1766;;; depending on whether x-i and y1 are the same.  We test this guess by
1767;;; determining whether guess*y2 is greater than the three high digits of x
1768;;; minus guess*y1 shifted left one digit:
1769;;;    ------------------------------
1770;;;   |    x-i    |   x-i-1  | x-i-2 |
1771;;;    ------------------------------
1772;;;    ------------------------------
1773;;; - | g*y1 high | g*y1 low |   0   |
1774;;;    ------------------------------
1775;;;                ...                   <   guess*y2     ???
1776;;; If guess*y2 is greater, then we decrement our guess by one and try again.
1777;;; This returns a guess that is either correct or one too large.
1778;;;
1779;;; the y's come from *truncate-y*, x's from *truncate-x*
1780;;; doing this in lap is not screamingly difficult - x's at i, i-1, i-2
1781
1782
1783
1784
1785
1786(defun bignum-truncate-guess-2 (x xidx y yidx)
1787  (digit-bind (guess-h guess-l)
1788              (%floor-99 x xidx y yidx)
1789    (truncate-guess-loop guess-h guess-l x xidx y yidx)))
1790
1791
1792
1793   
1794
1795;;; SHIFT-Y-FOR-TRUNCATE -- Internal.
1796;;;
1797;;; This returns the amount to shift y to place a one in the second highest
1798;;; bit.  Y must be positive.  If the last digit of y is zero, then y has a
1799;;; one in the previous digit's sign bit, so we know it will take one less
1800;;; than digit-size to get a one where we want.  Otherwise, we count how many
1801;;; right shifts it takes to get zero; subtracting this value from digit-size
1802;;; tells us how many high zeros there are which is one more than the shift
1803;;; amount sought.
1804;;;
1805;;; Note: This is exactly the same as one less than the integer-length of the
1806;;; last digit subtracted from the digit-size.
1807;;;
1808;;; We shift y to make it sufficiently large that doing the 64-bit by 32-bit
1809;;; %FLOOR calls ensures the quotient and remainder fit in 32-bits.
1810;;;
1811(defun shift-y-for-truncate (y)
1812  (the fixnum (1- (the fixnum (%bignum-sign-bits y)))))
1813
1814;;; SHIFT-AND-STORE-TRUNCATE-BUFFERS -- Internal.
1815;;;
1816;;; Stores two bignums into the truncation bignum buffers, shifting them on the
1817;;; way in.  This assumes x and y are positive and at least two in length, and
1818;;; it assumes *truncate-x* and *truncate-y* are one digit longer than x and y.
1819;;;
1820(defun shift-and-store-truncate-buffers (x len-x y len-y shift)
1821  (declare (type bignum-index len-x len-y)
1822           (type (integer 0 (#.digit-size)) shift))
1823  (cond ((eql 0 shift)
1824         (bignum-replace *truncate-x* x :end1 len-x)
1825         (bignum-replace *truncate-y* y :end1 len-y))
1826        (t
1827         (bignum-ashift-left-unaligned x 0 shift (the fixnum (1+ len-x)) *truncate-x*)
1828         (bignum-ashift-left-unaligned y 0 shift (the fixnum (1+ len-y)) *truncate-y*))))
1829
1830
1831
1832
1833;;;; General utilities.
1834
1835
1836;;; %NORMALIZE-BIGNUM-BUFFER -- Internal.
1837;;;
1838;;; Internal in-place operations use this to fixup remaining digits in the
1839;;; incoming data, such as in-place shifting.  This is basically the same as
1840;;; the first form in %NORMALIZE-BIGNUM, but we return the length of the buffer
1841;;; instead of shrinking the bignum.
1842;;;
1843
1844
1845
1846   
1847
1848
1849
1850
1851;;; %NORMALIZE-BIGNUM -- Internal.
1852;;;
1853;;; This drops the last digit if it is unnecessary sign information.  It
1854;;; repeats this as needed, possibly ending with a fixnum.  If the resulting
1855;;; length from shrinking is one, see if our one word is a fixnum.  Shift the
1856;;; possible fixnum bits completely out of the word, and compare this with
1857;;; shifting the sign bit all the way through.  If the bits are all 1's or 0's
1858;;; in both words, then there are just sign bits between the fixnum bits and
1859;;; the sign bit.  If we do have a fixnum, shift it over for the two low-tag
1860;;; bits.
1861;;;
1862
1863(defun %normalize-bignum (res)
1864  ;(declare (optimize (speed 3)(safety 0)))
1865  (%normalize-bignum-2 t res))
1866
1867;;; %MOSTLY-NORMALIZE-BIGNUM -- Internal.
1868;;;
1869;;; This drops the last digit if it is unnecessary sign information.  It
1870;;; repeats this as needed, possibly ending with a fixnum magnitude but never
1871;;; returning a fixnum.
1872;;;
1873
1874(defun %mostly-normalize-bignum (res &optional len)
1875  (declare (ignore len))
1876  (%normalize-bignum-2 nil res))
1877
1878
1879
1880
1881; think its ok
1882(defun ldb32 (hi-data lo-data size pos)
1883  (declare (fixnum hi-data lo-data size pos))
1884  (let* ((hi-bit (+ pos size))
1885         (mask (%i- (%ilsl size 1) 1)))
1886    (declare (fixnum hi-bit mask))   
1887    (%ilogand mask (if (< hi-bit 16)
1888                     (%iasr pos lo-data)
1889                     (if (>= pos 16)
1890                       (%ilsr (the fixnum (- pos 16)) hi-data)
1891                       (%ilogior 
1892                         (%iasr pos lo-data)
1893                         (%ilsl (the fixnum (- 16 pos)) hi-data)))))))
1894
1895
1896
1897
1898
1899; this was wrong for negative bigs when byte includes or exceeds sign
1900(defun %ldb-fixnum-from-bignum (bignum size position)
1901  (declare (fixnum size position))
1902  (let* ((low-idx (ash position -5))
1903         (low-bit (logand position 31))
1904         (hi-bit (+ low-bit size))
1905         (len (%bignum-length bignum))
1906         (minusp (bignum-minusp bignum)))
1907    (declare (fixnum size position low-bit hi-bit low-idx len))
1908    (if (>= low-idx len)
1909      (if minusp (1- (ash 1 size)) 0)     
1910      (multiple-value-bind (hi lo)(%bignum-ref bignum low-idx)
1911        (let ((chunk-lo (ldb32 hi lo (min size (%i- 32 low-bit)) low-bit)))
1912          (let ((val
1913                 (if (< hi-bit 32) 
1914                   chunk-lo
1915                   (progn
1916                     (setq low-idx (1+ low-idx))
1917                     (multiple-value-setq (hi lo)
1918                       (if (>= low-idx len)
1919                         (if minusp (values #xffff #xffff)(values 0 0))
1920                         (%bignum-ref bignum low-idx)))
1921                     (let ((chunk-hi (ldb32 hi lo (%i- size (%i- 32 low-bit)) 0)))
1922                       (%ilogior (ash chunk-hi (%i- 32 low-bit)) chunk-lo))))))
1923            val))))))
1924
1925(defun load-byte (size position integer)
1926  (if (and (bignump integer)
1927           (<= size (- 31 target::fixnumshift))
1928           (fixnump position))
1929    (%ldb-fixnum-from-bignum integer size position)
1930    (let ((mask (byte-mask size)))
1931      (if (and (fixnump mask) (fixnump integer)(fixnump position))
1932        (%ilogand mask (%iasr position integer))
1933        (logand mask (ash integer (- position)))))))   
1934
1935
1936#+safe-but-slow
1937(defun %bignum-bignum-gcd (u v)
1938  (setq u (abs u) v (abs v))
1939  (do* ((g 1 (ash g 1)))
1940       ((or (oddp u) (oddp v))
1941        (do* ()
1942             ((zerop u) (* g v))
1943          (cond ((evenp u) (setq u (ash u -1)))
1944                ((evenp v) (setq v (ash v -1)))
1945                (t (let* ((temp (ash (abs (- u v)) -1)))
1946                     (if (< u v)
1947                       (setq v temp)
1948                       (setq u temp)))))))
1949    (setq u (ash u -1) v (ash v -1))))
1950
1951
1952(defun %positive-bignum-bignum-gcd (u0 v0)
1953  (let* ((u-len (%bignum-length u0))
1954         (v-len (%bignum-length v0)))
1955    (declare (fixnum u-len v-len))
1956    (if (or (< u-len v-len)
1957            (and (= u-len v-len)
1958                 (< (bignum-compare u0 v0) 0)))
1959      (psetq u0 v0 v0 u0 u-len v-len v-len u-len))
1960    (with-bignum-buffers ((u u-len)
1961                          (u2 u-len)
1962                          (v v-len)
1963                          (v2 v-len))
1964      (bignum-replace u u0)
1965      (bignum-replace v v0)
1966      (let* ((u-trailing-0-bits (%bignum-count-trailing-zero-bits u))
1967             (u-trailing-0-digits (ash u-trailing-0-bits -5))
1968             (v-trailing-0-bits (%bignum-count-trailing-zero-bits v))
1969             (v-trailing-0-digits (ash v-trailing-0-bits -5)))
1970        (declare (fixnum u-trailing-0-bits v-trailing-0-bits))
1971        (unless (zerop u-trailing-0-bits)
1972          (bignum-shift-right-loop-1
1973           (logand u-trailing-0-bits 31)
1974           u2
1975           u
1976           (the fixnum (1- (the fixnum (- u-len u-trailing-0-digits ))))
1977           u-trailing-0-digits)
1978          (rotatef u u2)
1979          (%mostly-normalize-bignum-macro u)
1980          (setq u-len (%bignum-length u)))
1981        (unless (zerop v-trailing-0-bits)
1982          (bignum-shift-right-loop-1
1983           (logand v-trailing-0-bits 31)
1984           v2
1985           v
1986           (the fixnum (1- (the fixnum (- v-len v-trailing-0-digits))))
1987           v-trailing-0-digits)
1988          (rotatef v v2)
1989          (%mostly-normalize-bignum-macro v)
1990          (setq v-len (%bignum-length v)))
1991        (let* ((shift (min u-trailing-0-bits
1992                           v-trailing-0-bits)))
1993          (loop
1994            (let* ((fix-u (and (= u-len 1)
1995                               (let* ((hi-u (%bignum-ref-hi u 0)))
1996                                 (declare (fixnum hi-u))
1997                                 (= hi-u (the fixnum
1998                                           (logand hi-u (ash target::target-most-positive-fixnum -16)))))
1999                               (uvref u 0)))
2000                   (fix-v (and (= v-len 1)
2001                               (let* ((hi-v (%bignum-ref-hi v 0)))
2002                                 (declare (fixnum hi-v))
2003                                 (= hi-v (the fixnum
2004                                           (logand hi-v (ash target::target-most-positive-fixnum -16)))))
2005                               (uvref v 0))))
2006              (if fix-v
2007                (if fix-u
2008                  (return (ash (%fixnum-gcd fix-u fix-v) shift))
2009                  (return (ash (bignum-fixnum-gcd u fix-v) shift)))
2010                (if fix-u
2011                  (return (ash (bignum-fixnum-gcd v fix-u) shift)))))
2012             
2013            (let* ((signum (if (> u-len v-len)
2014                             1
2015                             (if (< u-len v-len)
2016                               -1
2017                               (bignum-compare u v)))))
2018              (declare (fixnum signum))
2019              (case signum
2020                (0                      ; (= u v)
2021                 (if (zerop shift)
2022                   (let* ((copy (%allocate-bignum u-len)))
2023                     (bignum-replace copy u)
2024                     (return copy))
2025                   (return (ash u shift))))
2026                (1                      ; (> u v)
2027                 (bignum-subtract-loop u u-len v v-len u)
2028                 (%mostly-normalize-bignum-macro u)
2029                 (setq u-len (%bignum-length u))
2030                 (setq u-trailing-0-bits
2031                       (%bignum-count-trailing-zero-bits u)
2032                       u-trailing-0-digits
2033                       (ash u-trailing-0-bits -5))
2034                 (unless (zerop u-trailing-0-bits)
2035                   (%init-misc 0 u2)
2036                   (bignum-shift-right-loop-1
2037                    (logand u-trailing-0-bits 31)
2038                    u2
2039                    u
2040                    (the fixnum (1- (the fixnum (- u-len
2041                                                   u-trailing-0-digits))))
2042                    u-trailing-0-digits)
2043                   (rotatef u u2)
2044                   (%mostly-normalize-bignum-macro u)
2045                   (setq u-len (%bignum-length u))))
2046                (t                      ; (> v u)
2047                 (bignum-subtract-loop v v-len u u-len v)
2048                 (%mostly-normalize-bignum-macro v)
2049                 (setq v-len (%bignum-length v))
2050                 (setq v-trailing-0-bits
2051                       (%bignum-count-trailing-zero-bits v)
2052                       v-trailing-0-digits
2053                       (ash v-trailing-0-bits -5))
2054                 (unless (zerop v-trailing-0-bits)
2055                   (%init-misc 0 v2)
2056                   (bignum-shift-right-loop-1
2057                    (logand v-trailing-0-bits 31)
2058                    v2
2059                    v
2060                    (the fixnum (1- (the fixnum (- v-len v-trailing-0-digits))))
2061                    v-trailing-0-digits)
2062                   (rotatef v v2)
2063                   (%mostly-normalize-bignum-macro v)
2064                   (setq v-len (%bignum-length v))))))))))))
2065
2066(defun %bignum-bignum-gcd (u v)
2067  (with-negated-bignum-buffers u v %positive-bignum-bignum-gcd))
2068
2069(defun unsignedwide->integer (uwidep)
2070  (with-bignum-buffers ((b 3))
2071    (setf (uvref b 0) (%get-unsigned-long uwidep 4)
2072          (uvref b 1) (%get-unsigned-long uwidep 0))
2073    (let* ((n (%normalize-bignum b)))
2074      (if (typep n 'bignum)
2075        (copy-bignum n)
2076        n))))
2077
2078(defun one-bignum-factor-of-two (a) 
2079  (declare (type bignum-type a))
2080  (let ((len (%bignum-length a)))
2081    (declare (fixnum len))
2082    (dotimes (i len)
2083      (multiple-value-bind (a-h a-l) (%bignum-ref a i)
2084        (declare (fixnum a-h a-l))
2085        (unless (and (= a-h 0)(= a-l 0))
2086          (return (+ (%ilsl 5 i)
2087                     (let* ((j 0)
2088                            (a a-l))
2089                       (declare (fixnum a j))
2090                       (if (= a-l 0) (setq j 16 a a-h))
2091                       (dotimes (i 16)           
2092                         (if (oddp a)
2093                           (return (%i+ j i))
2094                           (setq a (%iasr 1 a))))))))))))
2095
2096(defun logbitp (index integer)
2097  "Predicate returns T if bit index of integer is a 1."
2098  (number-case index
2099    (fixnum
2100     (if (minusp (the fixnum index))(report-bad-arg index '(integer 0))))
2101    (bignum
2102     ;; assuming bignum cant have more than most-positive-fixnum bits
2103     ;; (2 expt 24 longs)
2104     (if (bignum-minusp index)(report-bad-arg index '(integer 0)))
2105     ;; should error if integer isn't
2106     (return-from logbitp (minusp (require-type integer 'integer)))))
2107  (number-case integer
2108    (fixnum
2109     (if (%i< index (- target::nbits-in-word target::fixnumshift))
2110       (%ilogbitp index integer)
2111       (minusp (the fixnum integer))))
2112    (bignum
2113     (let ((bidx (%iasr 5 index))
2114           (bbit (%ilogand index 31)))
2115       (declare (fixnum bidx bbit))
2116       (if (>= bidx (%bignum-length integer))
2117         (bignum-minusp integer)
2118         (multiple-value-bind (hi lo) (%bignum-ref integer bidx)
2119           (declare (fixnum hi lo))
2120           (if (> bbit 15)
2121             (%ilogbitp (%i- bbit 16) hi)
2122             (%ilogbitp bbit lo))))))))
2123
2124) ; #+32-bit-target
Note: See TracBrowser for help on using the repository browser.