source: branches/working-0711/ccl/level-0/l0-bignum64.lisp @ 12973

Last change on this file since 12973 was 12973, checked in by gz, 10 years ago

Faster bignum multiplication (r12839, r12847, r12850)

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