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