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