src/code/bignum.lisp

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