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