[3139] | 1 | ;;-*- Mode: Lisp; Package: CCL -*- |
---|
| 2 | ;;; |
---|
[13067] | 3 | ;;; Copyright (C) 2009 Clozure Associates |
---|
[3139] | 4 | ;;; Copyright (C) 1994-2001 Digitool, Inc |
---|
[13066] | 5 | ;;; This file is part of Clozure CL. |
---|
[3139] | 6 | ;;; |
---|
[13066] | 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 |
---|
[3139] | 9 | ;;; file "LICENSE". The LLGPL consists of a preamble and the LGPL, |
---|
[13066] | 10 | ;;; which is distributed with Clozure CL as the file "LGPL". Where these |
---|
[3139] | 11 | ;;; conflict, the preamble takes precedence. |
---|
| 12 | ;;; |
---|
[13066] | 13 | ;;; Clozure CL is referenced in the preamble as the "LIBRARY." |
---|
[3139] | 14 | ;;; |
---|
| 15 | ;;; The LLGPL is also available online at |
---|
| 16 | ;;; http://opensource.franz.com/preamble.html |
---|
| 17 | |
---|
| 18 | (in-package "CCL") |
---|
| 19 | |
---|
[3140] | 20 | |
---|
[3139] | 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 | |
---|
[3140] | 130 | #+32-bit-target |
---|
| 131 | (progn |
---|
[3139] | 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 | |
---|
[11109] | 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))) |
---|
[3139] | 164 | |
---|
| 165 | |
---|
[11109] | 166 | |
---|
[3139] | 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) |
---|
[10157] | 688 | (>= len-b 16) |
---|
[14119] | 689 | #+(or x8632-target arm-target) |
---|
[10157] | 690 | nil) |
---|
[3139] | 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)) |
---|
[15578] | 733 | (%multiply-and-add-fixnum-loop bignum-len bignum fixnum result) |
---|
[3139] | 734 | (when negate-res |
---|
| 735 | (negate-bignum-in-place result)) |
---|
| 736 | (%normalize-bignum-macro result )))) |
---|
[9882] | 737 | (declare (dynamic-extent #'do-it)) |
---|
[3139] | 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)) |
---|
[14119] | 1461 | (with-bignum-buffers ((truncate-x len-x+1) |
---|
| 1462 | (truncate-y (the fixnum (1+ len-y)))) |
---|
[3139] | 1463 | (let ((y-shift (shift-y-for-truncate y))) |
---|
[14119] | 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) |
---|
[3139] | 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)) |
---|
[14119] | 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 ))) |
---|
[3139] | 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)) |
---|
[14119] | 1507 | (with-bignum-buffers ((truncate-x len-x+1) |
---|
| 1508 | (truncate-y (the fixnum (1+ len-y)))) |
---|
[3139] | 1509 | (let ((y-shift (shift-y-for-truncate y))) |
---|
[14119] | 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) |
---|
[3139] | 1512 | (when (not (eql 0 y-shift)) |
---|
| 1513 | (let* ((res-len-1 (1- len-y))) |
---|
| 1514 | (declare (fixnum res-len-1)) |
---|
[14119] | 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))) |
---|
[3139] | 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))))))))))))))) |
---|
[9882] | 1520 | (declare (dynamic-extent #'do-it)) |
---|
[3139] | 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)) |
---|
[9882] | 1603 | (declare (type bignum-index len-x)) |
---|
[3139] | 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)) |
---|
[9882] | 1625 | (declare (type bignum-index len-x)) |
---|
[3139] | 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 | |
---|
[14119] | 1653 | (defun do-truncate (truncate-x truncate-y len-x len-y) |
---|
[3139] | 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) |
---|
[14119] | 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)) |
---|
[3139] | 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 | |
---|
[14119] | 1681 | (defun do-truncate-no-quo (truncate-x truncate-y len-x len-y) |
---|
[3139] | 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 |
---|
[14119] | 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) |
---|
[3139] | 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 | |
---|
[14119] | 1710 | (defun try-bignum-truncate-guess (truncate-x truncate-y guess-h guess-l len-y low-x-digit) |
---|
[3139] | 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) |
---|
[14119] | 1721 | (multiple-value-bind (y-h y-l) (%bignum-ref truncate-y j) |
---|
[3139] | 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) |
---|
[14119] | 1731 | (multiple-value-bind (tx-h tx-l) (%bignum-ref truncate-x i) |
---|
[3139] | 1732 | (multiple-value-bind (x-h x-l temp-borrow) |
---|
| 1733 | (%subtract-with-borrow-1 tx-h tx-l low-h low-l borrow) |
---|
[14119] | 1734 | (%bignum-set truncate-x i x-h x-l) |
---|
[3139] | 1735 | (setq borrow temp-borrow))))) |
---|
| 1736 | (incf i)) |
---|
[14119] | 1737 | (multiple-value-bind (tx-h tx-l) (%bignum-ref truncate-x i) |
---|
[3139] | 1738 | (multiple-value-bind (x-h x-l) |
---|
| 1739 | (%subtract-with-borrow-1 tx-h tx-l carry-digit-h carry-digit-l borrow) |
---|
[14119] | 1740 | (%bignum-set truncate-x i x-h x-l))) |
---|
[3139] | 1741 | ;; See if guess is off by one, adding one Y back in if necessary. |
---|
| 1742 | |
---|
| 1743 | |
---|
[14119] | 1744 | (cond ((%digit-0-or-plusp truncate-x i) |
---|
[3139] | 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 |
---|
[14119] | 1750 | (bignum-add-loop-+ low-x-digit truncate-x truncate-y len-y) |
---|
[3139] | 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 | ;;; |
---|
[14119] | 1818 | (defun shift-and-store-truncate-buffers (truncate-x truncate-y x len-x y len-y shift) |
---|
[3139] | 1819 | (declare (type bignum-index len-x len-y) |
---|
| 1820 | (type (integer 0 (#.digit-size)) shift)) |
---|
| 1821 | (cond ((eql 0 shift) |
---|
[14119] | 1822 | (bignum-replace truncate-x x :end1 len-x) |
---|
| 1823 | (bignum-replace truncate-y y :end1 len-y)) |
---|
[3139] | 1824 | (t |
---|
[14119] | 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)))) |
---|
[3139] | 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)) |
---|
[5591] | 1969 | (unless (zerop u-trailing-0-bits) |
---|
[3139] | 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 |
---|
[5591] | 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 |
---|
[10157] | 1996 | (logand hi-u (ash target::target-most-positive-fixnum -16))))) |
---|
[5591] | 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 |
---|
[10157] | 2002 | (logand hi-v (ash target::target-most-positive-fixnum -16))))) |
---|
[5591] | 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))))) |
---|
[3139] | 2010 | |
---|
[5591] | 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) |
---|
[3139] | 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) |
---|
[5591] | 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) |
---|
[3139] | 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) |
---|
[5591] | 2062 | (setq v-len (%bignum-length v)))))))))))) |
---|
[3139] | 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 | |
---|
[9882] | 2122 | ) ; #+32-bit-target |
---|