source: trunk/source/level-0/l0-bignum64.lisp @ 8529

Last change on this file since 8529 was 8529, checked in by gb, 13 years ago

A couple of fixes in MULTIPLY-BIGNUM-AND-FIXNUM (including punting when
fixnum is most-negative-fixnum). Only broken on trunk,

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