src/code/bignum.lisp

   1 ;;;; code to implement bignum support
   2
   3 ;;;; This software is part of the SBCL system. See the README file for
   4 ;;;; more information.
   5 ;;;;
   6 ;;;; This software is derived from the CMU CL system, which was
   7 ;;;; written at Carnegie Mellon University and released into the
   8 ;;;; public domain. The software is in the public domain and is
   9 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
  10 ;;;; files for more information.
  11
  12 (in-package "SB!BIGNUM")
  13 \f
  14 ;;;; notes
  15
  16 ;;; comments from CMU CL:
  17 ;;;   These symbols define the interface to the number code:
  18 ;;;       add-bignums multiply-bignums negate-bignum subtract-bignum
  19 ;;;       multiply-bignum-and-fixnum multiply-fixnums
  20 ;;;       bignum-ashift-right bignum-ashift-left bignum-gcd
  21 ;;;       bignum-to-float bignum-integer-length
  22 ;;;       bignum-logical-and bignum-logical-ior bignum-logical-xor
  23 ;;;       bignum-logical-not bignum-load-byte bignum-deposit-byte
  24 ;;;       bignum-truncate bignum-plus-p bignum-compare make-small-bignum
  25 ;;;       bignum-logbitp bignum-logcount
  26 ;;;   These symbols define the interface to the compiler:
  27 ;;;       bignum-type bignum-element-type bignum-index %allocate-bignum
  28 ;;;       %bignum-length %bignum-set-length %bignum-ref %bignum-set
  29 ;;;       %digit-0-or-plusp %add-with-carry %subtract-with-borrow
  30 ;;;       %multiply-and-add %multiply %lognot %logand %logior %logxor
  31 ;;;       %fixnum-to-digit %floor %fixnum-digit-with-correct-sign %ashl
  32 ;;;       %ashr %digit-logical-shift-right))
  33
  34 ;;; The following interfaces will either be assembler routines or code
  35 ;;; sequences expanded into the code as basic bignum operations:
  36 ;;;    General:
  37 ;;;       %BIGNUM-LENGTH
  38 ;;;       %ALLOCATE-BIGNUM
  39 ;;;       %BIGNUM-REF
  40 ;;;       %NORMALIZE-BIGNUM
  41 ;;;       %BIGNUM-SET-LENGTH
  42 ;;;       %FIXNUM-DIGIT-WITH-CORRECT-SIGN
  43 ;;;       %SIGN-DIGIT
  44 ;;;       %ASHR
  45 ;;;       %ASHL
  46 ;;;       %BIGNUM-0-OR-PLUSP
  47 ;;;       %DIGIT-LOGICAL-SHIFT-RIGHT
  48 ;;;    General (May not exist when done due to sole use in %-routines.)
  49 ;;;       %DIGIT-0-OR-PLUSP
  50 ;;;    Addition:
  51 ;;;       %ADD-WITH-CARRY
  52 ;;;    Subtraction:
  53 ;;;       %SUBTRACT-WITH-BORROW
  54 ;;;    Multiplication
  55 ;;;       %MULTIPLY
  56 ;;;    Negation
  57 ;;;       %LOGNOT
  58 ;;;    Shifting (in place)
  59 ;;;       %NORMALIZE-BIGNUM-BUFFER
  60 ;;;    GCD/Relational operators:
  61 ;;;       %DIGIT-COMPARE
  62 ;;;       %DIGIT-GREATER
  63 ;;;    Relational operators:
  64 ;;;       %LOGAND
  65 ;;;       %LOGIOR
  66 ;;;       %LOGXOR
  67 ;;;    LDB
  68 ;;;       %FIXNUM-TO-DIGIT
  69 ;;;    TRUNCATE
  70 ;;;       %FLOOR
  71 ;;;
  72 ;;; Note: The floating routines know about the float representation.
  73 ;;;
  74 ;;; PROBLEM 1:
  75 ;;; There might be a problem with various LET's and parameters that take a
  76 ;;; digit value. We need to write these so those things stay in machine
  77 ;;; registers and number stack slots. I bind locals to these values, and I
  78 ;;; use function on them -- ZEROP, ASH, etc.
  79 ;;;
  80 ;;; PROBLEM 2:
  81 ;;; In shifting and byte operations, I use masks and logical operations that
  82 ;;; could result in intermediate bignums. This is hidden by the current system,
  83 ;;; but I may need to write these in a way that keeps these masks and logical
  84 ;;; operations from diving into the Lisp level bignum code.
  85 ;;;
  86 ;;; To do:
  87 ;;;    fixnums
  88 ;;;       logior, logxor, logand
  89 ;;;       depending on relationals, < (twice) and <= (twice)
  90 ;;;       or write compare thing (twice).
  91 ;;;       LDB on fixnum with bignum result.
  92 ;;;       DPB on fixnum with bignum result.
  93 ;;;       TRUNCATE returns zero or one as one value and fixnum or minus fixnum
  94 ;;;       for the other value when given (truncate fixnum bignum).
  95 ;;;       Returns (truncate bignum fixnum) otherwise.
  96 ;;;       addition
  97 ;;;       subtraction (twice)
  98 ;;;       multiply
  99 ;;;       GCD
 100 ;;;    Write MASK-FIELD and DEPOSIT-FIELD in terms of logical operations.
 101 ;;;    DIVIDE
 102 ;;;       IF (/ x y) with bignums:
 103 ;;;       do the truncate, and if rem is 0, return quotient.
 104 ;;;       if rem is non-0
 105 ;;;          gcd of x and y.
 106 ;;;          "truncate" each by gcd, ignoring remainder 0.
 107 ;;;          form ratio of each result, bottom is positive.
 108 \f
 109 ;;;; What's a bignum?
 110
 111 (defconstant digit-size sb!vm:n-word-bits)
 112
 113 (defconstant maximum-bignum-length (1- (ash 1 (- sb!vm:n-word-bits
 114                                                  sb!vm:n-widetag-bits))))
 115
 116 (defconstant all-ones-digit (1- (ash 1 sb!vm:n-word-bits)))
 117 \f
 118 ;;;; internal inline routines
 119
 120 ;;; %ALLOCATE-BIGNUM must zero all elements.
 121 (defun %allocate-bignum (length)
 122   (declare (type bignum-index length))
 123   (%allocate-bignum length))
 124
 125 ;;; Extract the length of the bignum.
 126 (defun %bignum-length (bignum)
 127   (declare (type bignum-type bignum))
 128   (%bignum-length bignum))
 129
 130 ;;; %BIGNUM-REF needs to access bignums as obviously as possible, and it needs
 131 ;;; to be able to return the digit somewhere no one looks for real objects.
 132 (defun %bignum-ref (bignum i)
 133   (declare (type bignum-type bignum)
 134            (type bignum-index i))
 135   (%bignum-ref bignum i))
 136 (defun %bignum-set (bignum i value)
 137   (declare (type bignum-type bignum)
 138            (type bignum-index i)
 139            (type bignum-element-type value))
 140   (%bignum-set bignum i value))
 141
 142 ;;; Return T if digit is positive, or NIL if negative.
 143 (defun %digit-0-or-plusp (digit)
 144   (declare (type bignum-element-type digit))
 145   (not (logbitp (1- digit-size) digit)))
 146
 147 #!-sb-fluid (declaim (inline %bignum-0-or-plusp))
 148 (defun %bignum-0-or-plusp (bignum len)
 149   (declare (type bignum-type bignum)
 150            (type bignum-index len))
 151   (%digit-0-or-plusp (%bignum-ref bignum (1- len))))
 152
 153 ;;; This should be in assembler, and should not cons intermediate
 154 ;;; results. It returns a bignum digit and a carry resulting from adding
 155 ;;; together a, b, and an incoming carry.
 156 (defun %add-with-carry (a b carry)
 157   (declare (type bignum-element-type a b)
 158            (type (mod 2) carry))
 159   (%add-with-carry a b carry))
 160
 161 ;;; This should be in assembler, and should not cons intermediate
 162 ;;; results. It returns a bignum digit and a borrow resulting from
 163 ;;; subtracting b from a, and subtracting a possible incoming borrow.
 164 ;;;
 165 ;;; We really do:  a - b - 1 + borrow, where borrow is either 0 or 1.
 166 (defun %subtract-with-borrow (a b borrow)
 167   (declare (type bignum-element-type a b)
 168            (type (mod 2) borrow))
 169   (%subtract-with-borrow a b borrow))
 170
 171 ;;; Multiply two digit-size numbers, returning a 2*digit-size result
 172 ;;; split into two digit-size quantities.
 173 (defun %multiply (x y)
 174   (declare (type bignum-element-type x y))
 175   (%multiply x y))
 176
 177 ;;; This multiplies x-digit and y-digit, producing high and low digits
 178 ;;; manifesting the result. Then it adds the low digit, res-digit, and
 179 ;;; carry-in-digit. Any carries (note, you still have to add two digits
 180 ;;; at a time possibly producing two carries) from adding these three
 181 ;;; digits get added to the high digit from the multiply, producing the
 182 ;;; next carry digit.  Res-digit is optional since two uses of this
 183 ;;; primitive multiplies a single digit bignum by a multiple digit
 184 ;;; bignum, and in this situation there is no need for a result buffer
 185 ;;; accumulating partial results which is where the res-digit comes
 186 ;;; from.
 187 (defun %multiply-and-add (x-digit y-digit carry-in-digit
 188                           &optional (res-digit 0))
 189   (declare (type bignum-element-type x-digit y-digit res-digit carry-in-digit))
 190   (%multiply-and-add x-digit y-digit carry-in-digit res-digit))
 191
 192 (defun %lognot (digit)
 193   (declare (type bignum-element-type digit))
 194   (%lognot digit))
 195
 196 ;;; Each of these does the digit-size unsigned op.
 197 #!-sb-fluid (declaim (inline %logand %logior %logxor))
 198 (defun %logand (a b)
 199   (declare (type bignum-element-type a b))
 200   (logand a b))
 201 (defun %logior (a b)
 202   (declare (type bignum-element-type a b))
 203   (logior a b))
 204 (defun %logxor (a b)
 205   (declare (type bignum-element-type a b))
 206   (logxor a b))
 207
 208 ;;; This takes a fixnum and sets it up as an unsigned digit-size
 209 ;;; quantity.
 210 (defun %fixnum-to-digit (x)
 211   (declare (fixnum x))
 212   (logand x (1- (ash 1 digit-size))))
 213
 214 #!-32x16-divide
 215 ;;; This takes three digits and returns the FLOOR'ed result of
 216 ;;; dividing the first two as a 2*digit-size integer by the third.
 217 ;;;
 218 ;;; Do weird LET and SETQ stuff to bamboozle the compiler into allowing
 219 ;;; the %FLOOR transform to expand into pseudo-assembler for which the
 220 ;;; compiler can later correctly allocate registers.
 221 (defun %floor (a b c)
 222   (let ((a a) (b b) (c c))
 223     (declare (type bignum-element-type a b c))
 224     (setq a a b b c c)
 225     (%floor a b c)))
 226
 227 ;;; Convert the digit to a regular integer assuming that the digit is signed.
 228 (defun %fixnum-digit-with-correct-sign (digit)
 229   (declare (type bignum-element-type digit))
 230   (if (logbitp (1- digit-size) digit)
 231       (logior digit (ash -1 digit-size))
 232       digit))
 233
 234 ;;; Do an arithmetic shift right of data even though bignum-element-type is
 235 ;;; unsigned.
 236 (defun %ashr (data count)
 237   (declare (type bignum-element-type data)
 238            (type (mod #.sb!vm:n-word-bits) count))
 239   (%ashr data count))
 240
 241 ;;; This takes a digit-size quantity and shifts it to the left,
 242 ;;; returning a digit-size quantity.
 243 (defun %ashl (data count)
 244   (declare (type bignum-element-type data)
 245            (type (mod #.sb!vm:n-word-bits) count))
 246   (%ashl data count))
 247
 248 ;;; Do an unsigned (logical) right shift of a digit by Count.
 249 (defun %digit-logical-shift-right (data count)
 250   (declare (type bignum-element-type data)
 251            (type (mod #.sb!vm:n-word-bits) count))
 252   (%digit-logical-shift-right data count))
 253
 254 ;;; Change the length of bignum to be newlen. Newlen must be the same or
 255 ;;; smaller than the old length, and any elements beyond newlen must be zeroed.
 256 (defun %bignum-set-length (bignum newlen)
 257   (declare (type bignum-type bignum)
 258            (type bignum-index newlen))
 259   (%bignum-set-length bignum newlen))
 260
 261 ;;; This returns 0 or "-1" depending on whether the bignum is positive. This
 262 ;;; is suitable for infinite sign extension to complete additions,
 263 ;;; subtractions, negations, etc. This cannot return a -1 represented as
 264 ;;; a negative fixnum since it would then have to low zeros.
 265 #!-sb-fluid (declaim (inline %sign-digit))
 266 (defun %sign-digit (bignum len)
 267   (declare (type bignum-type bignum)
 268            (type bignum-index len))
 269   (%ashr (%bignum-ref bignum (1- len)) (1- digit-size)))
 270
 271 ;;; These take two digit-size quantities and compare or contrast them
 272 ;;; without wasting time with incorrect type checking.
 273 #!-sb-fluid (declaim (inline %digit-compare %digit-greater))
 274 (defun %digit-compare (x y)
 275   (= x y))
 276 (defun %digit-greater (x y)
 277   (> x y))
 278 \f
 279 (declaim (optimize (speed 3) (safety 0)))
 280 \f
 281 ;;;; addition
 282
 283 (defun add-bignums (a b)
 284   (declare (type bignum-type a b))
 285   (let ((len-a (%bignum-length a))
 286         (len-b (%bignum-length b)))
 287     (declare (type bignum-index len-a len-b))
 288     (multiple-value-bind (a len-a b len-b)
 289         (if (> len-a len-b)
 290             (values a len-a b len-b)
 291             (values b len-b a len-a))
 292       (declare (type bignum-type a b)
 293                (type bignum-index len-a len-b))
 294       (let* ((len-res (1+ len-a))
 295              (res (%allocate-bignum len-res))
 296              (carry 0))
 297         (declare (type bignum-index len-res)
 298                  (type bignum-type res)
 299                  (type (mod 2) carry))
 300         (dotimes (i len-b)
 301           (declare (type bignum-index i))
 302           (multiple-value-bind (v k)
 303               (%add-with-carry (%bignum-ref a i) (%bignum-ref b i) carry)
 304             (declare (type bignum-element-type v)
 305                      (type (mod 2) k))
 306             (setf (%bignum-ref res i) v)
 307             (setf carry k)))
 308         (if (/= len-a len-b)
 309             (finish-add a res carry (%sign-digit b len-b) len-b len-a)
 310             (setf (%bignum-ref res len-a)
 311                   (%add-with-carry (%sign-digit a len-a)
 312                                    (%sign-digit b len-b)
 313                                    carry)))
 314         (%normalize-bignum res len-res)))))
 315
 316 ;;; This takes the longer of two bignums and propagates the carry through its
 317 ;;; remaining high order digits.
 318 (defun finish-add (a res carry sign-digit-b start end)
 319   (declare (type bignum-type a res)
 320            (type (mod 2) carry)
 321            (type bignum-element-type sign-digit-b)
 322            (type bignum-index start end))
 323   (do ((i start (1+ i)))
 324       ((= i end)
 325        (setf (%bignum-ref res end)
 326              (%add-with-carry (%sign-digit a end) sign-digit-b carry)))
 327     (declare (type bignum-index i))
 328     (multiple-value-bind (v k)
 329         (%add-with-carry (%bignum-ref a i) sign-digit-b carry)
 330       (setf (%bignum-ref res i) v)
 331       (setf carry k)))
 332   (values))
 333 \f
 334 ;;;; subtraction
 335
 336 (eval-when (:compile-toplevel :execute)
 337
 338 ;;; This subtracts b from a plugging result into res. Return-fun is the
 339 ;;; function to call that fixes up the result returning any useful values, such
 340 ;;; as the result. This macro may evaluate its arguments more than once.
 341 (sb!xc:defmacro subtract-bignum-loop (a len-a b len-b res len-res return-fun)
 342   (let ((borrow (gensym))
 343         (a-digit (gensym))
 344         (a-sign (gensym))
 345         (b-digit (gensym))
 346         (b-sign (gensym))
 347         (i (gensym))
 348         (v (gensym))
 349         (k (gensym)))
 350     `(let* ((,borrow 1)
 351             (,a-sign (%sign-digit ,a ,len-a))
 352             (,b-sign (%sign-digit ,b ,len-b)))
 353        (declare (type bignum-element-type ,a-sign ,b-sign))
 354        (dotimes (,i ,len-res)
 355          (declare (type bignum-index ,i))
 356          (let ((,a-digit (if (< ,i ,len-a) (%bignum-ref ,a ,i) ,a-sign))
 357                (,b-digit (if (< ,i ,len-b) (%bignum-ref ,b ,i) ,b-sign)))
 358            (declare (type bignum-element-type ,a-digit ,b-digit))
 359            (multiple-value-bind (,v ,k)
 360                (%subtract-with-borrow ,a-digit ,b-digit ,borrow)
 361              (setf (%bignum-ref ,res ,i) ,v)
 362              (setf ,borrow ,k))))
 363        (,return-fun ,res ,len-res))))
 364
 365 ) ;EVAL-WHEN
 366
 367 (defun subtract-bignum (a b)
 368   (declare (type bignum-type a b))
 369   (let* ((len-a (%bignum-length a))
 370          (len-b (%bignum-length b))
 371          (len-res (1+ (max len-a len-b)))
 372          (res (%allocate-bignum len-res)))
 373     (declare (type bignum-index len-a len-b len-res)) ;Test len-res for bounds?
 374     (subtract-bignum-loop a len-a b len-b res len-res %normalize-bignum)))
 375
 376 ;;; Operations requiring a subtraction without the overhead of intermediate
 377 ;;; results, such as GCD, use this. It assumes Result is big enough for the
 378 ;;; result.
 379 (defun subtract-bignum-buffers-with-len (a len-a b len-b result len-res)
 380   (declare (type bignum-type a b result)
 381            (type bignum-index len-a len-b len-res))
 382   (subtract-bignum-loop a len-a b len-b result len-res
 383                         %normalize-bignum-buffer))
 384
 385 (defun subtract-bignum-buffers (a len-a b len-b result)
 386   (declare (type bignum-type a b result)
 387            (type bignum-index len-a len-b))
 388   (subtract-bignum-loop a len-a b len-b result (max len-a len-b)
 389                         %normalize-bignum-buffer))
 390 \f
 391 ;;;; multiplication
 392
 393 (defun multiply-bignums (a b)
 394   (declare (type bignum-type a b))
 395   (let* ((a-plusp (%bignum-0-or-plusp a (%bignum-length a)))
 396          (b-plusp (%bignum-0-or-plusp b (%bignum-length b)))
 397          (a (if a-plusp a (negate-bignum a)))
 398          (b (if b-plusp b (negate-bignum b)))
 399          (len-a (%bignum-length a))
 400          (len-b (%bignum-length b))
 401          (len-res (+ len-a len-b))
 402          (res (%allocate-bignum len-res))
 403          (negate-res (not (eq a-plusp b-plusp))))
 404     (declare (type bignum-index len-a len-b len-res))
 405     (dotimes (i len-a)
 406       (declare (type bignum-index i))
 407       (let ((carry-digit 0)
 408             (x (%bignum-ref a i))
 409             (k i))
 410         (declare (type bignum-index k)
 411                  (type bignum-element-type carry-digit x))
 412         (dotimes (j len-b)
 413           (multiple-value-bind (big-carry res-digit)
 414               (%multiply-and-add x
 415                                  (%bignum-ref b j)
 416                                  (%bignum-ref res k)
 417                                  carry-digit)
 418             (declare (type bignum-element-type big-carry res-digit))
 419             (setf (%bignum-ref res k) res-digit)
 420             (setf carry-digit big-carry)
 421             (incf k)))
 422         (setf (%bignum-ref res k) carry-digit)))
 423     (when negate-res (negate-bignum-in-place res))
 424     (%normalize-bignum res len-res)))
 425
 426 (defun multiply-bignum-and-fixnum (bignum fixnum)
 427   (declare (type bignum-type bignum) (type fixnum fixnum))
 428   (let* ((bignum-plus-p (%bignum-0-or-plusp bignum (%bignum-length bignum)))
 429          (fixnum-plus-p (not (minusp fixnum)))
 430          (bignum (if bignum-plus-p bignum (negate-bignum bignum)))
 431          (bignum-len (%bignum-length bignum))
 432          (fixnum (if fixnum-plus-p fixnum (- fixnum)))
 433          (result (%allocate-bignum (1+ bignum-len)))
 434          (carry-digit 0))
 435     (declare (type bignum-type bignum result)
 436              (type bignum-index bignum-len)
 437              (type bignum-element-type fixnum carry-digit))
 438     (dotimes (index bignum-len)
 439       (declare (type bignum-index index))
 440       (multiple-value-bind (next-digit low)
 441           (%multiply-and-add (%bignum-ref bignum index) fixnum carry-digit)
 442         (declare (type bignum-element-type next-digit low))
 443         (setf carry-digit next-digit)
 444         (setf (%bignum-ref result index) low)))
 445     (setf (%bignum-ref result bignum-len) carry-digit)
 446     (unless (eq bignum-plus-p fixnum-plus-p)
 447       (negate-bignum-in-place result))
 448     (%normalize-bignum result (1+ bignum-len))))
 449
 450 (defun multiply-fixnums (a b)
 451   (declare (fixnum a b))
 452   (let* ((a-minusp (minusp a))
 453          (b-minusp (minusp b)))
 454     (multiple-value-bind (high low)
 455         (%multiply (if a-minusp (- a) a)
 456                    (if b-minusp (- b) b))
 457       (declare (type bignum-element-type high low))
 458       (if (and (zerop high)
 459                (%digit-0-or-plusp low))
 460           (let ((low (sb!ext:truly-the (unsigned-byte #.(1- sb!vm:n-word-bits))
 461                                        (%fixnum-digit-with-correct-sign low))))
 462             (if (eq a-minusp b-minusp)
 463                 low
 464                 (- low)))
 465           (let ((res (%allocate-bignum 2)))
 466             (%bignum-set res 0 low)
 467             (%bignum-set res 1 high)
 468             (unless (eq a-minusp b-minusp) (negate-bignum-in-place res))
 469             (%normalize-bignum res 2))))))
 470 \f
 471 ;;;; BIGNUM-REPLACE and WITH-BIGNUM-BUFFERS
 472
 473 (eval-when (:compile-toplevel :execute)
 474
 475 (sb!xc:defmacro bignum-replace (dest
 476                                 src
 477                                 &key
 478                                 (start1 '0)
 479                                 end1
 480                                 (start2 '0)
 481                                 end2
 482                                 from-end)
 483   (sb!int:once-only ((n-dest dest)
 484                      (n-src src))
 485     (let ((n-start1 (gensym))
 486           (n-end1 (gensym))
 487           (n-start2 (gensym))
 488           (n-end2 (gensym))
 489           (i1 (gensym))
 490           (i2 (gensym))
 491           (end1 (or end1 `(%bignum-length ,n-dest)))
 492           (end2 (or end2 `(%bignum-length ,n-src))))
 493       (if from-end
 494           `(let ((,n-start1 ,start1)
 495                  (,n-start2 ,start2))
 496              (do ((,i1 (1- ,end1) (1- ,i1))
 497                   (,i2 (1- ,end2) (1- ,i2)))
 498                  ((or (< ,i1 ,n-start1) (< ,i2 ,n-start2)))
 499                (declare (fixnum ,i1 ,i2))
 500                (%bignum-set ,n-dest ,i1
 501                             (%bignum-ref ,n-src ,i2))))
 502           `(let ((,n-end1 ,end1)
 503                  (,n-end2 ,end2))
 504              (do ((,i1 ,start1 (1+ ,i1))
 505                   (,i2 ,start2 (1+ ,i2)))
 506                  ((or (>= ,i1 ,n-end1) (>= ,i2 ,n-end2)))
 507                (declare (type bignum-index ,i1 ,i2))
 508                (%bignum-set ,n-dest ,i1
 509                             (%bignum-ref ,n-src ,i2))))))))
 510
 511 (sb!xc:defmacro with-bignum-buffers (specs &body body)
 512   #!+sb-doc
 513   "WITH-BIGNUM-BUFFERS ({(var size [init])}*) Form*"
 514   (sb!int:collect ((binds)
 515                    (inits))
 516     (dolist (spec specs)
 517       (let ((name (first spec))
 518             (size (second spec)))
 519         (binds `(,name (%allocate-bignum ,size)))
 520         (let ((init (third spec)))
 521           (when init
 522             (inits `(bignum-replace ,name ,init))))))
 523     `(let* ,(binds)
 524        ,@(inits)
 525        ,@body)))
 526
 527 ) ;EVAL-WHEN
 528 \f
 529 ;;;; GCD
 530
 531 (eval-when (:compile-toplevel :load-toplevel :execute)
 532   ;; The asserts in the GCD implementation are way too expensive to
 533   ;; check in normal use, and are disabled here.
 534   (sb!xc:defmacro gcd-assert (&rest args)
 535     (if nil
 536         `(assert ,@args)))
 537   ;; We'll be doing a lot of modular arithmetic.
 538   (sb!xc:defmacro modularly (form)
 539     `(logand all-ones-digit ,form)))
 540
 541 ;;; I'm not sure why I need this FTYPE declaration.  Compiled by the
 542 ;;; target compiler, it can deduce the return type fine, but without
 543 ;;; it, we pay a heavy price in BIGNUM-GCD when compiled by the
 544 ;;; cross-compiler. -- CSR, 2004-07-19
 545 (declaim (ftype (sfunction (bignum-type bignum-index bignum-type bignum-index)
 546                            sb!vm::positive-fixnum)
 547                 bignum-factors-of-two))
 548 (defun bignum-factors-of-two (a len-a b len-b)
 549   (declare (type bignum-index len-a len-b) (type bignum-type a b))
 550   (do ((i 0 (1+ i))
 551        (end (min len-a len-b)))
 552       ((= i end) (error "Unexpected zero bignums?"))
 553     (declare (type bignum-index i end))
 554     (let ((or-digits (%logior (%bignum-ref a i) (%bignum-ref b i))))
 555       (unless (zerop or-digits)
 556         (return (do ((j 0 (1+ j))
 557                      (or-digits or-digits (%ashr or-digits 1)))
 558                     ((oddp or-digits) (+ (* i digit-size) j))
 559                   (declare (type (mod #.sb!vm:n-word-bits) j))))))))
 560
 561 ;;; Multiply a bignum buffer with a fixnum or a digit, storing the
 562 ;;; result in another bignum buffer, and without using any
 563 ;;; temporaries. Inlined to avoid boxing smallnum if it's actually a
 564 ;;; digit. Needed by GCD, should possibly OAOO with
 565 ;;; MULTIPLY-BIGNUM-AND-FIXNUM.
 566 (declaim (inline multiply-bignum-buffer-and-smallnum-to-buffer))
 567 (defun multiply-bignum-buffer-and-smallnum-to-buffer (bignum bignum-len
 568                                                              smallnum res)
 569   (declare (type bignum-type bignum))
 570   (let* ((bignum-plus-p (%bignum-0-or-plusp bignum bignum-len))
 571          (smallnum-plus-p (not (minusp smallnum)))
 572          (smallnum (if smallnum-plus-p smallnum (- smallnum)))
 573          (carry-digit 0))
 574     (declare (type bignum-type bignum res)
 575              (type bignum-index bignum-len)
 576              (type bignum-element-type smallnum carry-digit))
 577     (unless bignum-plus-p
 578       (negate-bignum-buffer-in-place bignum bignum-len))
 579     (dotimes (index bignum-len)
 580       (declare (type bignum-index index))
 581       (multiple-value-bind (next-digit low)
 582           (%multiply-and-add (%bignum-ref bignum index)
 583                              smallnum
 584                              carry-digit)
 585         (declare (type bignum-element-type next-digit low))
 586         (setf carry-digit next-digit)
 587         (setf (%bignum-ref res index) low)))
 588     (setf (%bignum-ref res bignum-len) carry-digit)
 589     (unless bignum-plus-p
 590       (negate-bignum-buffer-in-place bignum bignum-len))
 591     (let ((res-len (%normalize-bignum-buffer res (1+ bignum-len))))
 592       (unless (eq bignum-plus-p smallnum-plus-p)
 593         (negate-bignum-buffer-in-place res res-len))
 594       res-len)))
 595
 596 ;;; Given U and V, return U / V mod 2^32. Implements the algorithm in the
 597 ;;; paper, but uses some clever bit-twiddling nicked from Nickle to do it.
 598 (declaim (inline bmod))
 599 (defun bmod (u v)
 600   (let ((ud (%bignum-ref u 0))
 601         (vd (%bignum-ref v 0))
 602         (umask 0)
 603         (imask 1)
 604         (m 0))
 605     (declare (type (unsigned-byte #.sb!vm:n-word-bits) ud vd umask imask m))
 606     (dotimes (i digit-size)
 607       (setf umask (logior umask imask))
 608       (when (logtest ud umask)
 609         (setf ud (modularly (- ud vd)))
 610         (setf m (modularly (logior m imask))))
 611       (setf imask (modularly (ash imask 1)))
 612       (setf vd (modularly (ash vd 1))))
 613     m))
 614
 615 (defun dmod (u u-len v v-len tmp1)
 616   (loop while (> (bignum-buffer-integer-length u u-len)
 617                  (+ (bignum-buffer-integer-length v v-len)
 618                     digit-size))
 619     do
 620     (unless (zerop (%bignum-ref u 0))
 621       (let* ((bmod (bmod u v))
 622              (tmp1-len (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
 623                                                                       bmod
 624                                                                       tmp1)))
 625         (setf u-len (subtract-bignum-buffers u u-len
 626                                              tmp1 tmp1-len
 627                                              u))
 628         (bignum-abs-buffer u u-len)))
 629     (gcd-assert (zerop (%bignum-ref u 0)))
 630     (setf u-len (bignum-buffer-ashift-right u u-len digit-size)))
 631   (let* ((d (+ 1 (- (bignum-buffer-integer-length u u-len)
 632                     (bignum-buffer-integer-length v v-len))))
 633          (n (1- (ash 1 d))))
 634     (declare (type (unsigned-byte #.(integer-length #.sb!vm:n-word-bits)) d)
 635              (type (unsigned-byte #.sb!vm:n-word-bits) n))
 636     (gcd-assert (>= d 0))
 637     (when (logtest (%bignum-ref u 0) n)
 638       (let ((tmp1-len
 639              (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
 640                                                             (logand n (bmod u
 641                                                                             v))
 642                                                             tmp1)))
 643         (setf u-len (subtract-bignum-buffers u u-len
 644                                              tmp1 tmp1-len
 645                                              u))
 646         (bignum-abs-buffer u u-len)))
 647     u-len))
 648
 649 (defconstant lower-ones-digit (1- (ash 1 (truncate sb!vm:n-word-bits 2))))
 650
 651 ;;; Find D and N such that (LOGAND ALL-ONES-DIGIT (- (* D X) (* N Y))) is 0,
 652 ;;; (< 0 N LOWER-ONES-DIGIT) and (< 0 (ABS D) LOWER-ONES-DIGIT).
 653 (defun reduced-ratio-mod (x y)
 654   (let* ((c (bmod x y))
 655          (n1 c)
 656          (d1 1)
 657          (n2 (modularly (1+ (modularly (lognot n1)))))
 658          (d2 (modularly -1)))
 659     (declare (type (unsigned-byte #.sb!vm:n-word-bits) n1 d1 n2 d2))
 660     (loop while (> n2 (expt 2 (truncate digit-size 2))) do
 661           (loop for i of-type (mod #.sb!vm:n-word-bits)
 662                 downfrom (- (integer-length n1) (integer-length n2))
 663                 while (>= n1 n2) do
 664                 (when (>= n1 (modularly (ash n2 i)))
 665                   (psetf n1 (modularly (- n1 (modularly (ash n2 i))))
 666                          d1 (modularly (- d1 (modularly (ash d2 i)))))))
 667           (psetf n1 n2
 668                  d1 d2
 669                  n2 n1
 670                  d2 d1))
 671     (values n2 (if (>= d2 (expt 2 (1- digit-size)))
 672                    (lognot (logand most-positive-fixnum (lognot d2)))
 673                    (logand lower-ones-digit d2)))))
 674
 675
 676 (defun copy-bignum (a &optional (len (%bignum-length a)))
 677   (let ((b (%allocate-bignum len)))
 678     (bignum-replace b a)
 679     (%bignum-set-length b len)
 680     b))
 681
 682 ;;; Allocate a single word bignum that holds fixnum. This is useful when
 683 ;;; we are trying to mix fixnum and bignum operands.
 684 #!-sb-fluid (declaim (inline make-small-bignum))
 685 (defun make-small-bignum (fixnum)
 686   (let ((res (%allocate-bignum 1)))
 687     (setf (%bignum-ref res 0) (%fixnum-to-digit fixnum))
 688     res))
 689
 690 ;; When the larger number is less than this many bignum digits long, revert
 691 ;; to old algorithm.
 692 (defparameter *accelerated-gcd-cutoff* 3)
 693
 694 ;;; Alternate between k-ary reduction with the help of
 695 ;;; REDUCED-RATIO-MOD and digit modulus reduction via DMOD. Once the
 696 ;;; arguments get small enough, drop through to BIGNUM-MOD-GCD (since
 697 ;;; k-ary reduction can introduce spurious factors, which need to be
 698 ;;; filtered out). Reference: Kenneth Weber, "The accelerated integer
 699 ;;; GCD algorithm", ACM Transactions on Mathematical Software, volume
 700 ;;; 21, number 1, March 1995, epp. 111-122.
 701 (defun bignum-gcd (u0 v0)
 702   (declare (type bignum-type u0 v0))
 703   (let* ((u1 (if (%bignum-0-or-plusp u0 (%bignum-length u0))
 704                  u0
 705                  (negate-bignum u0 nil)))
 706          (v1 (if (%bignum-0-or-plusp v0 (%bignum-length v0))
 707                  v0
 708                  (negate-bignum v0 nil))))
 709     (if (zerop v1)
 710         (return-from bignum-gcd u1))
 711     (when (> u1 v1)
 712       (rotatef u1 v1))
 713     (let ((n (mod v1 u1)))
 714       (setf v1 (if (fixnump n)
 715                    (make-small-bignum n)
 716                    n)))
 717     (if (and (= 1 (%bignum-length v1))
 718              (zerop (%bignum-ref v1 0)))
 719         (return-from bignum-gcd (%normalize-bignum u1
 720                                                    (%bignum-length u1))))
 721     (let* ((buffer-len (+ 2 (%bignum-length u1)))
 722            (u (%allocate-bignum buffer-len))
 723            (u-len (%bignum-length u1))
 724            (v (%allocate-bignum buffer-len))
 725            (v-len (%bignum-length v1))
 726            (tmp1 (%allocate-bignum buffer-len))
 727            (tmp1-len 0)
 728            (tmp2 (%allocate-bignum buffer-len))
 729            (tmp2-len 0)
 730            (factors-of-two
 731             (bignum-factors-of-two u1 (%bignum-length u1)
 732                                    v1 (%bignum-length v1))))
 733       (declare (type (or null bignum-index)
 734                      buffer-len u-len v-len tmp1-len tmp2-len))
 735       (bignum-replace u u1)
 736       (bignum-replace v v1)
 737       (setf u-len
 738             (make-gcd-bignum-odd u
 739                                  (bignum-buffer-ashift-right u u-len
 740                                                              factors-of-two)))
 741       (setf v-len
 742             (make-gcd-bignum-odd v
 743                                  (bignum-buffer-ashift-right v v-len
 744                                                              factors-of-two)))
 745       (loop until (or (< u-len *accelerated-gcd-cutoff*)
 746                       (not v-len)
 747                       (zerop v-len)
 748                       (and (= 1 v-len)
 749                            (zerop (%bignum-ref v 0))))
 750         do
 751         (gcd-assert (= buffer-len (%bignum-length u)
 752                        (%bignum-length v)
 753                        (%bignum-length tmp1)
 754                        (%bignum-length tmp2)))
 755         (if (> (bignum-buffer-integer-length u u-len)
 756                (+ #.(truncate sb!vm:n-word-bits 4)
 757                   (bignum-buffer-integer-length v v-len)))
 758             (setf u-len (dmod u u-len
 759                               v v-len
 760                               tmp1))
 761             (multiple-value-bind (n d) (reduced-ratio-mod u v)
 762               (setf tmp1-len
 763                     (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
 764                                                                    n tmp1))
 765               (setf tmp2-len
 766                     (multiply-bignum-buffer-and-smallnum-to-buffer u u-len
 767                                                                    d tmp2))
 768               (gcd-assert (= (copy-bignum tmp2 tmp2-len)
 769                              (* (copy-bignum u u-len) d)))
 770               (gcd-assert (= (copy-bignum tmp1 tmp1-len)
 771                              (* (copy-bignum v v-len) n)))
 772               (setf u-len
 773                     (subtract-bignum-buffers-with-len tmp1 tmp1-len
 774                                                       tmp2 tmp2-len
 775                                                       u
 776                                                       (1+ (max tmp1-len
 777                                                                tmp2-len))))
 778               (gcd-assert (or (zerop (- (copy-bignum tmp1 tmp1-len)
 779                                         (copy-bignum tmp2 tmp2-len)))
 780                               (= (copy-bignum u u-len)
 781                                  (- (copy-bignum tmp1 tmp1-len)
 782                                     (copy-bignum tmp2 tmp2-len)))))
 783               (bignum-abs-buffer u u-len)
 784               (gcd-assert (zerop (modularly u)))))
 785         (setf u-len (make-gcd-bignum-odd u u-len))
 786         (rotatef u v)
 787         (rotatef u-len v-len))
 788       (setf u (copy-bignum u u-len))
 789       (let ((n (bignum-mod-gcd v1 u)))
 790         (ash (bignum-mod-gcd u1 (if (fixnump n)
 791                                     (make-small-bignum n)
 792                                     n))
 793              factors-of-two)))))
 794
 795 (defun bignum-mod-gcd (a b)
 796   (declare (type bignum-type a b))
 797   (when (< a b)
 798     (rotatef a b))
 799   ;; While the length difference of A and B is sufficiently large,
 800   ;; reduce using MOD (slowish, but it should equalize the sizes of
 801   ;; A and B pretty quickly). After that, use the binary GCD
 802   ;; algorithm to handle the rest.
 803   (loop until (and (= (%bignum-length b) 1) (zerop (%bignum-ref b 0))) do
 804         (when (<= (%bignum-length a) (1+ (%bignum-length b)))
 805           (return-from bignum-mod-gcd (bignum-binary-gcd a b)))
 806         (let ((rem (mod a b)))
 807           (if (fixnump rem)
 808               (setf a (make-small-bignum rem))
 809               (setf a rem))
 810           (rotatef a b)))
 811   (if (= (%bignum-length a) 1)
 812       (%normalize-bignum a 1)
 813       a))
 814
 815 (defun bignum-binary-gcd (a b)
 816   (declare (type bignum-type a b))
 817   (let* ((len-a (%bignum-length a))
 818          (len-b (%bignum-length b)))
 819     (declare (type bignum-index len-a len-b))
 820     (with-bignum-buffers ((a-buffer len-a a)
 821                           (b-buffer len-b b)
 822                           (res-buffer (max len-a len-b)))
 823       (let* ((factors-of-two
 824               (bignum-factors-of-two a-buffer len-a
 825                                      b-buffer len-b))
 826              (len-a (make-gcd-bignum-odd
 827                      a-buffer
 828                      (bignum-buffer-ashift-right a-buffer len-a
 829                                                  factors-of-two)))
 830              (len-b (make-gcd-bignum-odd
 831                      b-buffer
 832                      (bignum-buffer-ashift-right b-buffer len-b
 833                                                  factors-of-two))))
 834         (declare (type bignum-index len-a len-b))
 835         (let ((x a-buffer)
 836               (len-x len-a)
 837               (y b-buffer)
 838               (len-y len-b)
 839               (z res-buffer))
 840           (loop
 841             (multiple-value-bind (u v len-v r len-r)
 842                 (bignum-gcd-order-and-subtract x len-x y len-y z)
 843               (declare (type bignum-index len-v len-r))
 844               (when (and (= len-r 1) (zerop (%bignum-ref r 0)))
 845                 (if (zerop factors-of-two)
 846                     (let ((ret (%allocate-bignum len-v)))
 847                       (dotimes (i len-v)
 848                         (setf (%bignum-ref ret i) (%bignum-ref v i)))
 849                       (return (%normalize-bignum ret len-v)))
 850                     (return (bignum-ashift-left v factors-of-two len-v))))
 851               (setf x v  len-x len-v)
 852               (setf y r  len-y (make-gcd-bignum-odd r len-r))
 853               (setf z u))))))))
 854
 855 (defun bignum-gcd-order-and-subtract (a len-a b len-b res)
 856   (declare (type bignum-index len-a len-b) (type bignum-type a b))
 857   (cond ((= len-a len-b)
 858          (do ((i (1- len-a) (1- i)))
 859              ((= i -1)
 860               (setf (%bignum-ref res 0) 0)
 861               (values a b len-b res 1))
 862            (let ((a-digit (%bignum-ref a i))
 863                  (b-digit (%bignum-ref b i)))
 864              (cond ((%digit-compare a-digit b-digit))
 865                    ((%digit-greater a-digit b-digit)
 866                     (return
 867                      (values a b len-b res
 868                              (subtract-bignum-buffers a len-a b len-b
 869                                                       res))))
 870                    (t
 871                     (return
 872                      (values b a len-a res
 873                              (subtract-bignum-buffers b len-b
 874                                                       a len-a
 875                                                       res))))))))
 876         ((> len-a len-b)
 877          (values a b len-b res
 878                  (subtract-bignum-buffers a len-a b len-b res)))
 879         (t
 880          (values b a len-a res
 881                  (subtract-bignum-buffers b len-b a len-a res)))))
 882
 883 (defun make-gcd-bignum-odd (a len-a)
 884   (declare (type bignum-type a) (type bignum-index len-a))
 885   (dotimes (index len-a)
 886     (declare (type bignum-index index))
 887     (do ((digit (%bignum-ref a index) (%ashr digit 1))
 888          (increment 0 (1+ increment)))
 889         ((zerop digit))
 890       (declare (type (mod #.sb!vm:n-word-bits) increment))
 891       (when (oddp digit)
 892         (return-from make-gcd-bignum-odd
 893                      (bignum-buffer-ashift-right a len-a
 894                                                  (+ (* index digit-size)
 895                                                     increment)))))))
 896
 897 \f
 898 ;;;; negation
 899
 900 (eval-when (:compile-toplevel :execute)
 901
 902 ;;; This negates bignum-len digits of bignum, storing the resulting digits into
 903 ;;; result (possibly EQ to bignum) and returning whatever end-carry there is.
 904 (sb!xc:defmacro bignum-negate-loop (bignum
 905                                     bignum-len
 906                                     &optional (result nil resultp))
 907   (let ((carry (gensym))
 908         (end (gensym))
 909         (value (gensym))
 910         (last (gensym)))
 911     `(let* (,@(if (not resultp) `(,last))
 912             (,carry
 913              (multiple-value-bind (,value ,carry)
 914                  (%add-with-carry (%lognot (%bignum-ref ,bignum 0)) 1 0)
 915                ,(if resultp
 916                     `(setf (%bignum-ref ,result 0) ,value)
 917                     `(setf ,last ,value))
 918                ,carry))
 919             (i 1)
 920             (,end ,bignum-len))
 921        (declare (type bit ,carry)
 922                 (type bignum-index i ,end))
 923        (loop
 924          (when (= i ,end) (return))
 925          (multiple-value-bind (,value temp)
 926              (%add-with-carry (%lognot (%bignum-ref ,bignum i)) 0 ,carry)
 927            ,(if resultp
 928                 `(setf (%bignum-ref ,result i) ,value)
 929                 `(setf ,last ,value))
 930            (setf ,carry temp))
 931          (incf i))
 932        ,(if resultp carry `(values ,carry ,last)))))
 933
 934 ) ; EVAL-WHEN
 935
 936 ;;; Fully-normalize is an internal optional. It cause this to always return
 937 ;;; a bignum, without any extraneous digits, and it never returns a fixnum.
 938 (defun negate-bignum (x &optional (fully-normalize t))
 939   (declare (type bignum-type x))
 940   (let* ((len-x (%bignum-length x))
 941          (len-res (1+ len-x))
 942          (res (%allocate-bignum len-res)))
 943     (declare (type bignum-index len-x len-res)) ;Test len-res for range?
 944     (let ((carry (bignum-negate-loop x len-x res)))
 945       (setf (%bignum-ref res len-x)
 946             (%add-with-carry (%lognot (%sign-digit x len-x)) 0 carry)))
 947     (if fully-normalize
 948         (%normalize-bignum res len-res)
 949         (%mostly-normalize-bignum res len-res))))
 950
 951 ;;; This assumes bignum is positive; that is, the result of negating it will
 952 ;;; stay in the provided allocated bignum.
 953 (defun negate-bignum-buffer-in-place (bignum bignum-len)
 954   (bignum-negate-loop bignum bignum-len bignum)
 955   bignum)
 956
 957 (defun negate-bignum-in-place (bignum)
 958   (declare (inline negate-bignum-buffer-in-place))
 959   (negate-bignum-buffer-in-place bignum (%bignum-length bignum)))
 960
 961 (defun bignum-abs-buffer (bignum len)
 962   (unless (%bignum-0-or-plusp bignum len)
 963     (negate-bignum-buffer-in-place bignum len)))
 964 \f
 965 ;;;; shifting
 966
 967 (eval-when (:compile-toplevel :execute)
 968
 969 ;;; This macro is used by BIGNUM-ASHIFT-RIGHT, BIGNUM-BUFFER-ASHIFT-RIGHT, and
 970 ;;; BIGNUM-LDB-BIGNUM-RES. They supply a termination form that references
 971 ;;; locals established by this form. Source is the source bignum. Start-digit
 972 ;;; is the first digit in source from which we pull bits. Start-pos is the
 973 ;;; first bit we want. Res-len-form is the form that computes the length of
 974 ;;; the resulting bignum. Termination is a DO termination form with a test and
 975 ;;; body. When result is supplied, it is the variable to which this binds a
 976 ;;; newly allocated bignum.
 977 ;;;
 978 ;;; Given start-pos, 1-31 inclusively, of shift, we form the j'th resulting
 979 ;;; digit from high bits of the i'th source digit and the start-pos number of
 980 ;;; bits from the i+1'th source digit.
 981 (sb!xc:defmacro shift-right-unaligned (source
 982                                        start-digit
 983                                        start-pos
 984                                        res-len-form
 985                                        termination
 986                                        &optional result)
 987   `(let* ((high-bits-in-first-digit (- digit-size ,start-pos))
 988           (res-len ,res-len-form)
 989           (res-len-1 (1- res-len))
 990           ,@(if result `((,result (%allocate-bignum res-len)))))
 991      (declare (type bignum-index res-len res-len-1))
 992      (do ((i ,start-digit i+1)
 993           (i+1 (1+ ,start-digit) (1+ i+1))
 994           (j 0 (1+ j)))
 995          ,termination
 996        (declare (type bignum-index i i+1 j))
 997        (setf (%bignum-ref ,(if result result source) j)
 998              (%logior (%digit-logical-shift-right (%bignum-ref ,source i)
 999                                                   ,start-pos)
1000                       (%ashl (%bignum-ref ,source i+1)
1001                              high-bits-in-first-digit))))))
1002
1003 ) ; EVAL-WHEN
1004
1005 ;;; First compute the number of whole digits to shift, shifting them by
1006 ;;; skipping them when we start to pick up bits, and the number of bits to
1007 ;;; shift the remaining digits into place. If the number of digits is greater
1008 ;;; than the length of the bignum, then the result is either 0 or -1. If we
1009 ;;; shift on a digit boundary (that is, n-bits is zero), then we just copy
1010 ;;; digits. The last branch handles the general case which uses a macro that a
1011 ;;; couple other routines use. The fifth argument to the macro references
1012 ;;; locals established by the macro.
1013 (defun bignum-ashift-right (bignum count)
1014   (declare (type bignum-type bignum)
1015            (type unsigned-byte count))
1016   (let ((bignum-len (%bignum-length bignum)))
1017     (declare (type bignum-index bignum-len))
1018     (cond ((fixnump count)
1019            (multiple-value-bind (digits n-bits) (truncate count digit-size)
1020              (declare (type bignum-index digits))
1021              (cond
1022               ((>= digits bignum-len)
1023                (if (%bignum-0-or-plusp bignum bignum-len) 0 -1))
1024               ((zerop n-bits)
1025                (bignum-ashift-right-digits bignum digits))
1026               (t
1027                (shift-right-unaligned bignum digits n-bits (- bignum-len digits)
1028                                       ((= j res-len-1)
1029                                        (setf (%bignum-ref res j)
1030                                              (%ashr (%bignum-ref bignum i) n-bits))
1031                                        (%normalize-bignum res res-len))
1032                                       res)))))
1033           ((> count bignum-len)
1034            (if (%bignum-0-or-plusp bignum bignum-len) 0 -1))
1035            ;; Since a FIXNUM should be big enough to address anything in
1036            ;; memory, including arrays of bits, and since arrays of bits
1037            ;; take up about the same space as corresponding fixnums, there
1038            ;; should be no way that we fall through to this case: any shift
1039            ;; right by a bignum should give zero. But let's check anyway:
1040           (t (error "bignum overflow: can't shift right by ~S" count)))))
1041
1042 (defun bignum-ashift-right-digits (bignum digits)
1043   (declare (type bignum-type bignum)
1044            (type bignum-index digits))
1045   (let* ((res-len (- (%bignum-length bignum) digits))
1046          (res (%allocate-bignum res-len)))
1047     (declare (type bignum-index res-len)
1048              (type bignum-type res))
1049     (bignum-replace res bignum :start2 digits)
1050     (%normalize-bignum res res-len)))
1051
1052 ;;; GCD uses this for an in-place shifting operation. This is different enough
1053 ;;; from BIGNUM-ASHIFT-RIGHT that it isn't worth folding the bodies into a
1054 ;;; macro, but they share the basic algorithm. This routine foregoes a first
1055 ;;; test for digits being greater than or equal to bignum-len since that will
1056 ;;; never happen for its uses in GCD. We did fold the last branch into a macro
1057 ;;; since it was duplicated a few times, and the fifth argument to it
1058 ;;; references locals established by the macro.
1059 (defun bignum-buffer-ashift-right (bignum bignum-len x)
1060   (declare (type bignum-index bignum-len) (fixnum x))
1061   (multiple-value-bind (digits n-bits) (truncate x digit-size)
1062     (declare (type bignum-index digits))
1063     (cond
1064      ((zerop n-bits)
1065       (let ((new-end (- bignum-len digits)))
1066         (bignum-replace bignum bignum :end1 new-end :start2 digits
1067                         :end2 bignum-len)
1068         (%normalize-bignum-buffer bignum new-end)))
1069      (t
1070       (shift-right-unaligned bignum digits n-bits (- bignum-len digits)
1071                              ((= j res-len-1)
1072                               (setf (%bignum-ref bignum j)
1073                                     (%ashr (%bignum-ref bignum i) n-bits))
1074                               (%normalize-bignum-buffer bignum res-len)))))))
1075
1076 ;;; This handles shifting a bignum buffer to provide fresh bignum data for some
1077 ;;; internal routines. We know bignum is safe when called with bignum-len.
1078 ;;; First we compute the number of whole digits to shift, shifting them
1079 ;;; starting to store farther along the result bignum. If we shift on a digit
1080 ;;; boundary (that is, n-bits is zero), then we just copy digits. The last
1081 ;;; branch handles the general case.
1082 (defun bignum-ashift-left (bignum x &optional bignum-len)
1083   (declare (type bignum-type bignum)
1084            (type unsigned-byte x)
1085            (type (or null bignum-index) bignum-len))
1086   (if (fixnump x)
1087     (multiple-value-bind (digits n-bits) (truncate x digit-size)
1088       (let* ((bignum-len (or bignum-len (%bignum-length bignum)))
1089              (res-len (+ digits bignum-len 1)))
1090         (when (> res-len maximum-bignum-length)
1091           (error "can't represent result of left shift"))
1092         (if (zerop n-bits)
1093           (bignum-ashift-left-digits bignum bignum-len digits)
1094           (bignum-ashift-left-unaligned bignum digits n-bits res-len))))
1095     ;; Left shift by a number too big to be represented as a fixnum
1096     ;; would exceed our memory capacity, since a fixnum is big enough
1097     ;; to index any array, including a bit array.
1098     (error "can't represent result of left shift")))
1099
1100 (defun bignum-ashift-left-digits (bignum bignum-len digits)
1101   (declare (type bignum-index bignum-len digits))
1102   (let* ((res-len (+ bignum-len digits))
1103          (res (%allocate-bignum res-len)))
1104     (declare (type bignum-index res-len))
1105     (bignum-replace res bignum :start1 digits :end1 res-len :end2 bignum-len
1106                     :from-end t)
1107     res))
1108
1109 ;;; BIGNUM-TRUNCATE uses this to store into a bignum buffer by supplying res.
1110 ;;; When res comes in non-nil, then this foregoes allocating a result, and it
1111 ;;; normalizes the buffer instead of the would-be allocated result.
1112 ;;;
1113 ;;; We start storing into one digit higher than digits, storing a whole result
1114 ;;; digit from parts of two contiguous digits from bignum. When the loop
1115 ;;; finishes, we store the remaining bits from bignum's first digit in the
1116 ;;; first non-zero result digit, digits. We also grab some left over high
1117 ;;; bits from the last digit of bignum.
1118 (defun bignum-ashift-left-unaligned (bignum digits n-bits res-len
1119                                      &optional (res nil resp))
1120   (declare (type bignum-index digits res-len)
1121            (type (mod #.digit-size) n-bits))
1122   (let* ((remaining-bits (- digit-size n-bits))
1123          (res-len-1 (1- res-len))
1124          (res (or res (%allocate-bignum res-len))))
1125     (declare (type bignum-index res-len res-len-1))
1126     (do ((i 0 i+1)
1127          (i+1 1 (1+ i+1))
1128          (j (1+ digits) (1+ j)))
1129         ((= j res-len-1)
1130          (setf (%bignum-ref res digits)
1131                (%ashl (%bignum-ref bignum 0) n-bits))
1132          (setf (%bignum-ref res j)
1133                (%ashr (%bignum-ref bignum i) remaining-bits))
1134          (if resp
1135              (%normalize-bignum-buffer res res-len)
1136              (%normalize-bignum res res-len)))
1137       (declare (type bignum-index i i+1 j))
1138       (setf (%bignum-ref res j)
1139             (%logior (%digit-logical-shift-right (%bignum-ref bignum i)
1140                                                  remaining-bits)
1141                      (%ashl (%bignum-ref bignum i+1) n-bits))))))
1142 \f
1143 ;;;; relational operators
1144
1145 ;;; Return T iff bignum is positive.
1146 (defun bignum-plus-p (bignum)
1147   (declare (type bignum-type bignum))
1148   (%bignum-0-or-plusp bignum (%bignum-length bignum)))
1149
1150 ;;; This compares two bignums returning -1, 0, or 1, depending on
1151 ;;; whether a is less than, equal to, or greater than b.
1152 (declaim (ftype (function (bignum bignum) (integer -1 1)) bignum-compare))
1153 (defun bignum-compare (a b)
1154   (declare (type bignum-type a b))
1155   (let* ((len-a (%bignum-length a))
1156          (len-b (%bignum-length b))
1157          (a-plusp (%bignum-0-or-plusp a len-a))
1158          (b-plusp (%bignum-0-or-plusp b len-b)))
1159     (declare (type bignum-index len-a len-b))
1160     (cond ((not (eq a-plusp b-plusp))
1161            (if a-plusp 1 -1))
1162           ((= len-a len-b)
1163            (do ((i (1- len-a) (1- i)))
1164                (())
1165              (declare (type bignum-index i))
1166              (let ((a-digit (%bignum-ref a i))
1167                    (b-digit (%bignum-ref b i)))
1168                (declare (type bignum-element-type a-digit b-digit))
1169                (when (%digit-greater a-digit b-digit)
1170                  (return 1))
1171                (when (%digit-greater b-digit a-digit)
1172                  (return -1)))
1173              (when (zerop i) (return 0))))
1174           ((> len-a len-b)
1175            (if a-plusp 1 -1))
1176           (t (if a-plusp -1 1)))))
1177 \f
1178 ;;;; float conversion
1179
1180 ;;; Make a single or double float with the specified significand,
1181 ;;; exponent and sign.
1182 (defun single-float-from-bits (bits exp plusp)
1183   (declare (fixnum exp))
1184   (declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
1185   (let ((res (dpb exp
1186                   sb!vm:single-float-exponent-byte
1187                   (logandc2 (logand #xffffffff
1188                                     (%bignum-ref bits 1))
1189                             sb!vm:single-float-hidden-bit))))
1190     (make-single-float
1191      (if plusp
1192          res
1193          (logior res (ash -1 sb!vm:float-sign-shift))))))
1194 (defun double-float-from-bits (bits exp plusp)
1195   (declare (fixnum exp))
1196   (declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
1197   (let ((hi (dpb exp
1198                  sb!vm:double-float-exponent-byte
1199                  (logandc2 (ecase sb!vm::n-word-bits
1200                              (32 (%bignum-ref bits 2))
1201                              (64 (ash (%bignum-ref bits 1) -32)))
1202                            sb!vm:double-float-hidden-bit)))
1203         (lo (logand #xffffffff (%bignum-ref bits 1))))
1204     (make-double-float (if plusp
1205                            hi
1206                            (logior hi (ash -1 sb!vm:float-sign-shift)))
1207                        lo)))
1208 #!+(and long-float x86)
1209 (defun long-float-from-bits (bits exp plusp)
1210   (declare (fixnum exp))
1211   (declare (optimize #-sb-xc-host (sb!ext:inhibit-warnings 3)))
1212   (make-long-float
1213    (if plusp
1214        exp
1215        (logior exp (ash 1 15)))
1216    (%bignum-ref bits 2)
1217    (%bignum-ref bits 1)))
1218
1219 ;;; Convert Bignum to a float in the specified Format, rounding to the best
1220 ;;; approximation.
1221 (defun bignum-to-float (bignum format)
1222   (let* ((plusp (bignum-plus-p bignum))
1223          (x (if plusp bignum (negate-bignum bignum)))
1224          (len (bignum-integer-length x))
1225          (digits (float-format-digits format))
1226          (keep (+ digits digit-size))
1227          (shift (- keep len))
1228          (shifted (if (minusp shift)
1229                       (bignum-ashift-right x (- shift))
1230                       (bignum-ashift-left x shift)))
1231          (low (%bignum-ref shifted 0))
1232          (round-bit (ash 1 (1- digit-size))))
1233     (declare (type bignum-index len digits keep) (fixnum shift))
1234     (labels ((round-up ()
1235                (let ((rounded (add-bignums shifted round-bit)))
1236                  (if (> (integer-length rounded) keep)
1237                      (float-from-bits (bignum-ashift-right rounded 1)
1238                                       (1+ len))
1239                      (float-from-bits rounded len))))
1240              (float-from-bits (bits len)
1241                (declare (type bignum-index len))
1242                (ecase format
1243                  (single-float
1244                   (single-float-from-bits
1245                    bits
1246                    (check-exponent len sb!vm:single-float-bias
1247                                    sb!vm:single-float-normal-exponent-max)
1248                    plusp))
1249                  (double-float
1250                   (double-float-from-bits
1251                    bits
1252                    (check-exponent len sb!vm:double-float-bias
1253                                    sb!vm:double-float-normal-exponent-max)
1254                    plusp))
1255                  #!+long-float
1256                  (long-float
1257                   (long-float-from-bits
1258                    bits
1259                    (check-exponent len sb!vm:long-float-bias
1260                                    sb!vm:long-float-normal-exponent-max)
1261                    plusp))))
1262              (check-exponent (exp bias max)
1263                (declare (type bignum-index len))
1264                (let ((exp (+ exp bias)))
1265                  (when (> exp max)
1266                    ;; Why a SIMPLE-TYPE-ERROR? Well, this is mainly
1267                    ;; called by COERCE, which requires an error of
1268                    ;; TYPE-ERROR if the conversion can't happen
1269                    ;; (except in certain circumstances when we are
1270                    ;; coercing to a FUNCTION) -- CSR, 2002-09-18
1271                    (error 'simple-type-error
1272                           :format-control "Too large to be represented as a ~S:~%  ~S"
1273                           :format-arguments (list format x)
1274                           :expected-type format
1275                           :datum x))
1276                  exp)))
1277
1278     (cond
1279      ;; Round down if round bit is 0.
1280      ((not (logtest round-bit low))
1281       (float-from-bits shifted len))
1282      ;; If only round bit is set, then round to even.
1283      ((and (= low round-bit)
1284            (dotimes (i (- (%bignum-length x) (ceiling keep digit-size))
1285                        t)
1286              (unless (zerop (%bignum-ref x i)) (return nil))))
1287       (let ((next (%bignum-ref shifted 1)))
1288         (if (oddp next)
1289             (round-up)
1290             (float-from-bits shifted len))))
1291      ;; Otherwise, round up.
1292      (t
1293       (round-up))))))
1294 \f
1295 ;;;; integer length and logbitp/logcount
1296
1297 (defun bignum-buffer-integer-length (bignum len)
1298   (declare (type bignum-type bignum))
1299   (let* ((len-1 (1- len))
1300          (digit (%bignum-ref bignum len-1)))
1301     (declare (type bignum-index len len-1)
1302              (type bignum-element-type digit))
1303     (+ (integer-length (%fixnum-digit-with-correct-sign digit))
1304        (* len-1 digit-size))))
1305
1306 (defun bignum-integer-length (bignum)
1307   (declare (type bignum-type bignum))
1308   (bignum-buffer-integer-length bignum (%bignum-length bignum)))
1309
1310 (defun bignum-logbitp (index bignum)
1311   (declare (type bignum-type bignum))
1312   (let ((len (%bignum-length bignum)))
1313     (declare (type bignum-index len))
1314     (multiple-value-bind (word-index bit-index)
1315         (floor index digit-size)
1316       (if (>= word-index len)
1317           (not (bignum-plus-p bignum))
1318           (logbitp bit-index (%bignum-ref bignum word-index))))))
1319
1320 (defun bignum-logcount (bignum)
1321   (declare (type bignum-type bignum))
1322   (let ((length (%bignum-length bignum))
1323         (result 0))
1324     (declare (type bignum-index length)
1325              (fixnum result))
1326     (do ((index 0 (1+ index)))
1327         ((= index length)
1328          (if (%bignum-0-or-plusp bignum length)
1329              result
1330              (- (* length digit-size) result)))
1331       (let ((digit (%bignum-ref bignum index)))
1332         (declare (type bignum-element-type digit))
1333         (incf result (logcount digit))))))
1334 \f
1335 ;;;; logical operations
1336
1337 ;;;; NOT
1338
1339 (defun bignum-logical-not (a)
1340   (declare (type bignum-type a))
1341   (let* ((len (%bignum-length a))
1342          (res (%allocate-bignum len)))
1343     (declare (type bignum-index len))
1344     (dotimes (i len res)
1345       (declare (type bignum-index i))
1346       (setf (%bignum-ref res i) (%lognot (%bignum-ref a i))))))
1347
1348 ;;;; AND
1349
1350 (defun bignum-logical-and (a b)
1351   (declare (type bignum-type a b))
1352   (let* ((len-a (%bignum-length a))
1353          (len-b (%bignum-length b))
1354          (a-plusp (%bignum-0-or-plusp a len-a))
1355          (b-plusp (%bignum-0-or-plusp b len-b)))
1356     (declare (type bignum-index len-a len-b))
1357     (cond
1358      ((< len-a len-b)
1359       (if a-plusp
1360           (logand-shorter-positive a len-a b (%allocate-bignum len-a))
1361           (logand-shorter-negative a len-a b len-b (%allocate-bignum len-b))))
1362      ((< len-b len-a)
1363       (if b-plusp
1364           (logand-shorter-positive b len-b a (%allocate-bignum len-b))
1365           (logand-shorter-negative b len-b a len-a (%allocate-bignum len-a))))
1366      (t (logand-shorter-positive a len-a b (%allocate-bignum len-a))))))
1367
1368 ;;; This takes a shorter bignum, a and len-a, that is positive. Because this
1369 ;;; is AND, we don't care about any bits longer than a's since its infinite 0
1370 ;;; sign bits will mask the other bits out of b. The result is len-a big.
1371 (defun logand-shorter-positive (a len-a b res)
1372   (declare (type bignum-type a b res)
1373            (type bignum-index len-a))
1374   (dotimes (i len-a)
1375     (declare (type bignum-index i))
1376     (setf (%bignum-ref res i)
1377           (%logand (%bignum-ref a i) (%bignum-ref b i))))
1378   (%normalize-bignum res len-a))
1379
1380 ;;; This takes a shorter bignum, a and len-a, that is negative. Because this
1381 ;;; is AND, we just copy any bits longer than a's since its infinite 1 sign
1382 ;;; bits will include any bits from b. The result is len-b big.
1383 (defun logand-shorter-negative (a len-a b len-b res)
1384   (declare (type bignum-type a b res)
1385            (type bignum-index len-a len-b))
1386   (dotimes (i len-a)
1387     (declare (type bignum-index i))
1388     (setf (%bignum-ref res i)
1389           (%logand (%bignum-ref a i) (%bignum-ref b i))))
1390   (do ((i len-a (1+ i)))
1391       ((= i len-b))
1392     (declare (type bignum-index i))
1393     (setf (%bignum-ref res i) (%bignum-ref b i)))
1394   (%normalize-bignum res len-b))
1395
1396 ;;;; IOR
1397
1398 (defun bignum-logical-ior (a b)
1399   (declare (type bignum-type a b))
1400   (let* ((len-a (%bignum-length a))
1401          (len-b (%bignum-length b))
1402          (a-plusp (%bignum-0-or-plusp a len-a))
1403          (b-plusp (%bignum-0-or-plusp b len-b)))
1404     (declare (type bignum-index len-a len-b))
1405     (cond
1406      ((< len-a len-b)
1407       (if a-plusp
1408           (logior-shorter-positive a len-a b len-b (%allocate-bignum len-b))
1409           (logior-shorter-negative a len-a b len-b (%allocate-bignum len-b))))
1410      ((< len-b len-a)
1411       (if b-plusp
1412           (logior-shorter-positive b len-b a len-a (%allocate-bignum len-a))
1413           (logior-shorter-negative b len-b a len-a (%allocate-bignum len-a))))
1414      (t (logior-shorter-positive a len-a b len-b (%allocate-bignum len-a))))))
1415
1416 ;;; This takes a shorter bignum, a and len-a, that is positive. Because this
1417 ;;; is IOR, we don't care about any bits longer than a's since its infinite
1418 ;;; 0 sign bits will mask the other bits out of b out to len-b. The result
1419 ;;; is len-b long.
1420 (defun logior-shorter-positive (a len-a b len-b res)
1421   (declare (type bignum-type a b res)
1422            (type bignum-index len-a len-b))
1423   (dotimes (i len-a)
1424     (declare (type bignum-index i))
1425     (setf (%bignum-ref res i)
1426           (%logior (%bignum-ref a i) (%bignum-ref b i))))
1427   (do ((i len-a (1+ i)))
1428       ((= i len-b))
1429     (declare (type bignum-index i))
1430     (setf (%bignum-ref res i) (%bignum-ref b i)))
1431   (%normalize-bignum res len-b))
1432
1433 ;;; This takes a shorter bignum, a and len-a, that is negative. Because this
1434 ;;; is IOR, we just copy any bits longer than a's since its infinite 1 sign
1435 ;;; bits will include any bits from b. The result is len-b long.
1436 (defun logior-shorter-negative (a len-a b len-b res)
1437   (declare (type bignum-type a b res)
1438            (type bignum-index len-a len-b))
1439   (dotimes (i len-a)
1440     (declare (type bignum-index i))
1441     (setf (%bignum-ref res i)
1442           (%logior (%bignum-ref a i) (%bignum-ref b i))))
1443   (do ((i len-a (1+ i))
1444        (sign (%sign-digit a len-a)))
1445       ((= i len-b))
1446     (declare (type bignum-index i))
1447     (setf (%bignum-ref res i) sign))
1448   (%normalize-bignum res len-b))
1449
1450 ;;;; XOR
1451
1452 (defun bignum-logical-xor (a b)
1453   (declare (type bignum-type a b))
1454   (let ((len-a (%bignum-length a))
1455         (len-b (%bignum-length b)))
1456     (declare (type bignum-index len-a len-b))
1457     (if (< len-a len-b)
1458         (bignum-logical-xor-aux a len-a b len-b (%allocate-bignum len-b))
1459         (bignum-logical-xor-aux b len-b a len-a (%allocate-bignum len-a)))))
1460
1461 ;;; This takes the shorter of two bignums in a and len-a. Res is len-b
1462 ;;; long. Do the XOR.
1463 (defun bignum-logical-xor-aux (a len-a b len-b res)
1464   (declare (type bignum-type a b res)
1465            (type bignum-index len-a len-b))
1466   (dotimes (i len-a)
1467     (declare (type bignum-index i))
1468     (setf (%bignum-ref res i)
1469           (%logxor (%bignum-ref a i) (%bignum-ref b i))))
1470   (do ((i len-a (1+ i))
1471        (sign (%sign-digit a len-a)))
1472       ((= i len-b))
1473     (declare (type bignum-index i))
1474     (setf (%bignum-ref res i) (%logxor sign (%bignum-ref b i))))
1475   (%normalize-bignum res len-b))
1476 \f
1477 ;;;; There used to be a bunch of code to implement "efficient" versions of LDB
1478 ;;;; and DPB here.  But it apparently was never used, so it's been deleted.
1479 ;;;;   --njf, 2007-02-04
1480 \f
1481 ;;;; TRUNCATE
1482
1483 ;;; This is the original sketch of the algorithm from which I implemented this
1484 ;;; TRUNCATE, assuming both operands are bignums. I should modify this to work
1485 ;;; with the documentation on my functions, as a general introduction. I've
1486 ;;; left this here just in case someone needs it in the future. Don't look at
1487 ;;; this unless reading the functions' comments leaves you at a loss. Remember
1488 ;;; this comes from Knuth, so the book might give you the right general
1489 ;;; overview.
1490 ;;;
1491 ;;; (truncate x y):
1492 ;;;
1493 ;;; If X's magnitude is less than Y's, then result is 0 with remainder X.
1494 ;;;
1495 ;;; Make x and y positive, copying x if it is already positive.
1496 ;;;
1497 ;;; Shift y left until there's a 1 in the 30'th bit (most significant, non-sign
1498 ;;;       digit)
1499 ;;;    Just do most sig digit to determine how much to shift whole number.
1500 ;;; Shift x this much too.
1501 ;;; Remember this initial shift count.
1502 ;;;
1503 ;;; Allocate q to be len-x minus len-y quantity plus 1.
1504 ;;;
1505 ;;; i = last digit of x.
1506 ;;; k = last digit of q.
1507 ;;;
1508 ;;; LOOP
1509 ;;;
1510 ;;; j = last digit of y.
1511 ;;;
1512 ;;; compute guess.
1513 ;;; if x[i] = y[j] then g = (1- (ash 1 digit-size))
1514 ;;; else g = x[i]x[i-1]/y[j].
1515 ;;;
1516 ;;; check guess.
1517 ;;; %UNSIGNED-MULTIPLY returns b and c defined below.
1518 ;;;    a = x[i-1] - (logand (* g y[j]) #xFFFFFFFF).
1519 ;;;       Use %UNSIGNED-MULTIPLY taking low-order result.
1520 ;;;    b = (logand (ash (* g y[j-1]) (- digit-size)) (1- (ash 1 digit-size))).
1521 ;;;    c = (logand (* g y[j-1]) (1- (ash 1 digit-size))).
1522 ;;; if a < b, okay.
1523 ;;; if a > b, guess is too high
1524 ;;;    g = g - 1; go back to "check guess".
1525 ;;; if a = b and c > x[i-2], guess is too high
1526 ;;;    g = g - 1; go back to "check guess".
1527 ;;; GUESS IS 32-BIT NUMBER, SO USE THING TO KEEP IN SPECIAL REGISTER
1528 ;;; SAME FOR A, B, AND C.
1529 ;;;
1530 ;;; Subtract g * y from x[i - len-y+1]..x[i]. See paper for doing this in step.
1531 ;;; If x[i] < 0, guess is screwed up.
1532 ;;;    negative g, then add 1
1533 ;;;    zero or positive g, then subtract 1
1534 ;;; AND add y back into x[len-y+1..i].
1535 ;;;
1536 ;;; q[k] = g.
1537 ;;; i = i - 1.
1538 ;;; k = k - 1.
1539 ;;;
1540 ;;; If k>=0, goto LOOP.
1541 ;;;
1542 ;;; Now quotient is good, but remainder is not.
1543 ;;; Shift x right by saved initial left shifting count.
1544 ;;;
1545 ;;; Check quotient and remainder signs.
1546 ;;; x pos y pos --> q pos r pos
1547 ;;; x pos y neg --> q neg r pos
1548 ;;; x neg y pos --> q neg r neg
1549 ;;; x neg y neg --> q pos r neg
1550 ;;;
1551 ;;; Normalize quotient and remainder. Cons result if necessary.
1552
1553
1554 ;;; This used to be split into multiple functions, which shared state
1555 ;;; in special variables *TRUNCATE-X* and *TRUNCATE-Y*. Having so many
1556 ;;; special variable accesses in tight inner loops was having a large
1557 ;;; effect on performance, so the helper functions have now been
1558 ;;; refactored into local functions and the special variables into
1559 ;;; lexicals.  There was also a lot of boxing and unboxing of
1560 ;;; (UNSIGNED-BYTE 32)'s going on, which this refactoring
1561 ;;; eliminated. This improves the performance on some CL-BENCH tests
1562 ;;; by up to 50%, which is probably signigicant enough to justify the
1563 ;;; reduction in readability that was introduced. --JES, 2004-08-07
1564 (defun bignum-truncate (x y)
1565   (declare (type bignum-type x y))
1566   (let (truncate-x truncate-y)
1567     (labels
1568         ;;; Divide X by Y when Y is a single bignum digit. BIGNUM-TRUNCATE
1569         ;;; fixes up the quotient and remainder with respect to sign and
1570         ;;; normalization.
1571         ;;;
1572         ;;; We don't have to worry about shifting Y to make its most
1573         ;;; significant digit sufficiently large for %FLOOR to return
1574         ;;; digit-size quantities for the q-digit and r-digit. If Y is
1575         ;;; a single digit bignum, it is already large enough for
1576         ;;; %FLOOR. That is, it has some bits on pretty high in the
1577         ;;; digit.
1578         ((bignum-truncate-single-digit (x len-x y)
1579            (declare (type bignum-index len-x))
1580            (let ((q (%allocate-bignum len-x))
1581                  (r 0)
1582                  (y (%bignum-ref y 0)))
1583              (declare (type bignum-element-type r y))
1584              (do ((i (1- len-x) (1- i)))
1585                  ((minusp i))
1586                (multiple-value-bind (q-digit r-digit)
1587                    (%floor r (%bignum-ref x i) y)
1588                  (declare (type bignum-element-type q-digit r-digit))
1589                  (setf (%bignum-ref q i) q-digit)
1590                  (setf r r-digit)))
1591              (let ((rem (%allocate-bignum 1)))
1592                (setf (%bignum-ref rem 0) r)
1593                (values q rem))))
1594         ;;; This returns a guess for the next division step. Y1 is the
1595         ;;; highest y digit, and y2 is the second to highest y
1596         ;;; digit. The x... variables are the three highest x digits
1597         ;;; for the next division step.
1598         ;;;
1599         ;;; From Knuth, our guess is either all ones or x-i and x-i-1
1600         ;;; divided by y1, depending on whether x-i and y1 are the
1601         ;;; same. We test this guess by determining whether guess*y2
1602         ;;; is greater than the three high digits of x minus guess*y1
1603         ;;; shifted left one digit:
1604         ;;;    ------------------------------
1605         ;;;   |    x-i    |   x-i-1  | x-i-2 |
1606         ;;;    ------------------------------
1607         ;;;    ------------------------------
1608         ;;; - | g*y1 high | g*y1 low |   0   |
1609         ;;;    ------------------------------
1610         ;;;             ...               <   guess*y2     ???
1611         ;;; If guess*y2 is greater, then we decrement our guess by one
1612         ;;; and try again.  This returns a guess that is either
1613         ;;; correct or one too large.
1614          (bignum-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
1615            (declare (type bignum-element-type y1 y2 x-i x-i-1 x-i-2))
1616            (let ((guess (if (%digit-compare x-i y1)
1617                             all-ones-digit
1618                             (%floor x-i x-i-1 y1))))
1619              (declare (type bignum-element-type guess))
1620              (loop
1621                  (multiple-value-bind (high-guess*y1 low-guess*y1)
1622                      (%multiply guess y1)
1623                    (declare (type bignum-element-type low-guess*y1
1624                                   high-guess*y1))
1625                    (multiple-value-bind (high-guess*y2 low-guess*y2)
1626                        (%multiply guess y2)
1627                      (declare (type bignum-element-type high-guess*y2
1628                                     low-guess*y2))
1629                      (multiple-value-bind (middle-digit borrow)
1630                          (%subtract-with-borrow x-i-1 low-guess*y1 1)
1631                        (declare (type bignum-element-type middle-digit)
1632                                 (fixnum borrow))
1633                        ;; Supplying borrow of 1 means there was no
1634                        ;; borrow, and we know x-i-2 minus 0 requires
1635                        ;; no borrow.
1636                        (let ((high-digit (%subtract-with-borrow x-i
1637                                                                 high-guess*y1
1638                                                                 borrow)))
1639                          (declare (type bignum-element-type high-digit))
1640                          (if (and (%digit-compare high-digit 0)
1641                                   (or (%digit-greater high-guess*y2
1642                                                       middle-digit)
1643                                       (and (%digit-compare middle-digit
1644                                                            high-guess*y2)
1645                                            (%digit-greater low-guess*y2
1646                                                            x-i-2))))
1647                              (setf guess (%subtract-with-borrow guess 1 1))
1648                              (return guess)))))))))
1649         ;;; Divide TRUNCATE-X by TRUNCATE-Y, returning the quotient
1650         ;;; and destructively modifying TRUNCATE-X so that it holds
1651         ;;; the remainder.
1652         ;;;
1653         ;;; LEN-X and LEN-Y tell us how much of the buffers we care about.
1654         ;;;
1655         ;;; TRUNCATE-X definitely has at least three digits, and it has one
1656         ;;; more than TRUNCATE-Y. This keeps i, i-1, i-2, and low-x-digit
1657         ;;; happy. Thanks to SHIFT-AND-STORE-TRUNCATE-BUFFERS.
1658          (return-quotient-leaving-remainder (len-x len-y)
1659            (declare (type bignum-index len-x len-y))
1660            (let* ((len-q (- len-x len-y))
1661                   ;; Add one for extra sign digit in case high bit is on.
1662                   (q (%allocate-bignum (1+ len-q)))
1663                   (k (1- len-q))
1664                   (y1 (%bignum-ref truncate-y (1- len-y)))
1665                   (y2 (%bignum-ref truncate-y (- len-y 2)))
1666                   (i (1- len-x))
1667                   (i-1 (1- i))
1668                   (i-2 (1- i-1))
1669                   (low-x-digit (- i len-y)))
1670              (declare (type bignum-index len-q k i i-1 i-2 low-x-digit)
1671                       (type bignum-element-type y1 y2))
1672              (loop
1673                  (setf (%bignum-ref q k)
1674                        (try-bignum-truncate-guess
1675                         ;; This modifies TRUNCATE-X. Must access
1676                         ;; elements each pass.
1677                         (bignum-truncate-guess y1 y2
1678                                                (%bignum-ref truncate-x i)
1679                                                (%bignum-ref truncate-x i-1)
1680                                                (%bignum-ref truncate-x i-2))
1681                         len-y low-x-digit))
1682                  (cond ((zerop k) (return))
1683                        (t (decf k)
1684                           (decf low-x-digit)
1685                           (shiftf i i-1 i-2 (1- i-2)))))
1686              q))
1687         ;;; This takes a digit guess, multiplies it by TRUNCATE-Y for a
1688         ;;; result one greater in length than LEN-Y, and subtracts this result
1689         ;;; from TRUNCATE-X. LOW-X-DIGIT is the first digit of X to start
1690         ;;; the subtraction, and we know X is long enough to subtract a LEN-Y
1691         ;;; plus one length bignum from it. Next we check the result of the
1692         ;;; subtraction, and if the high digit in X became negative, then our
1693         ;;; guess was one too big. In this case, return one less than GUESS
1694         ;;; passed in, and add one value of Y back into X to account for
1695         ;;; subtracting one too many. Knuth shows that the guess is wrong on
1696         ;;; the order of 3/b, where b is the base (2 to the digit-size power)
1697         ;;; -- pretty rarely.
1698          (try-bignum-truncate-guess (guess len-y low-x-digit)
1699            (declare (type bignum-index low-x-digit len-y)
1700                     (type bignum-element-type guess))
1701            (let ((carry-digit 0)
1702                  (borrow 1)
1703                  (i low-x-digit))
1704              (declare (type bignum-element-type carry-digit)
1705                       (type bignum-index i)
1706                       (fixnum borrow))
1707              ;; Multiply guess and divisor, subtracting from dividend
1708              ;; simultaneously.
1709              (dotimes (j len-y)
1710                (multiple-value-bind (high-digit low-digit)
1711                    (%multiply-and-add guess
1712                                       (%bignum-ref truncate-y j)
1713                                       carry-digit)
1714                  (declare (type bignum-element-type high-digit low-digit))
1715                  (setf carry-digit high-digit)
1716                  (multiple-value-bind (x temp-borrow)
1717                      (%subtract-with-borrow (%bignum-ref truncate-x i)
1718                                             low-digit
1719                                             borrow)
1720                    (declare (type bignum-element-type x)
1721                             (fixnum temp-borrow))
1722                    (setf (%bignum-ref truncate-x i) x)
1723                    (setf borrow temp-borrow)))
1724                (incf i))
1725              (setf (%bignum-ref truncate-x i)
1726                    (%subtract-with-borrow (%bignum-ref truncate-x i)
1727                                           carry-digit borrow))
1728              ;; See whether guess is off by one, adding one
1729              ;; Y back in if necessary.
1730              (cond ((%digit-0-or-plusp (%bignum-ref truncate-x i))
1731                     guess)
1732                    (t
1733                     ;; If subtraction has negative result, add one
1734                     ;; divisor value back in. The guess was one too
1735                     ;; large in magnitude.
1736                     (let ((i low-x-digit)
1737                           (carry 0))
1738                       (dotimes (j len-y)
1739                         (multiple-value-bind (v k)
1740                             (%add-with-carry (%bignum-ref truncate-y j)
1741                                              (%bignum-ref truncate-x i)
1742                                              carry)
1743                           (declare (type bignum-element-type v))
1744                           (setf (%bignum-ref truncate-x i) v)
1745                           (setf carry k))
1746                         (incf i))
1747                       (setf (%bignum-ref truncate-x i)
1748                             (%add-with-carry (%bignum-ref truncate-x i)
1749                                              0 carry)))
1750                     (%subtract-with-borrow guess 1 1)))))
1751         ;;; This returns the amount to shift y to place a one in the
1752         ;;; second highest bit. Y must be positive. If the last digit
1753         ;;; of y is zero, then y has a one in the previous digit's
1754         ;;; sign bit, so we know it will take one less than digit-size
1755         ;;; to get a one where we want. Otherwise, we count how many
1756         ;;; right shifts it takes to get zero; subtracting this value
1757         ;;; from digit-size tells us how many high zeros there are
1758         ;;; which is one more than the shift amount sought.
1759         ;;;
1760         ;;; Note: This is exactly the same as one less than the
1761         ;;; integer-length of the last digit subtracted from the
1762         ;;; digit-size.
1763         ;;;
1764         ;;; We shift y to make it sufficiently large that doing the
1765         ;;; 2*digit-size by digit-size %FLOOR calls ensures the quotient and
1766         ;;; remainder fit in digit-size.
1767          (shift-y-for-truncate (y)
1768            (let* ((len (%bignum-length y))
1769                   (last (%bignum-ref y (1- len))))
1770              (declare (type bignum-index len)
1771                       (type bignum-element-type last))
1772              (- digit-size (integer-length last) 1)))
1773          ;;; Stores two bignums into the truncation bignum buffers,
1774          ;;; shifting them on the way in. This assumes x and y are
1775          ;;; positive and at least two in length, and it assumes
1776          ;;; truncate-x and truncate-y are one digit longer than x and
1777          ;;; y.
1778          (shift-and-store-truncate-buffers (x len-x y len-y shift)
1779            (declare (type bignum-index len-x len-y)
1780                     (type (integer 0 (#.digit-size)) shift))
1781            (cond ((zerop shift)
1782                   (bignum-replace truncate-x x :end1 len-x)
1783                   (bignum-replace truncate-y y :end1 len-y))
1784                  (t
1785                   (bignum-ashift-left-unaligned x 0 shift (1+ len-x)
1786                                                 truncate-x)
1787                   (bignum-ashift-left-unaligned y 0 shift (1+ len-y)
1788                                                 truncate-y))))) ;; LABELS
1789       ;;; Divide X by Y returning the quotient and remainder. In the
1790       ;;; general case, we shift Y to set up for the algorithm, and we
1791       ;;; use two buffers to save consing intermediate values. X gets
1792       ;;; destructively modified to become the remainder, and we have
1793       ;;; to shift it to account for the initial Y shift. After we
1794       ;;; multiple bind q and r, we first fix up the signs and then
1795       ;;; return the normalized results.
1796       (let* ((x-plusp (%bignum-0-or-plusp x (%bignum-length x)))
1797              (y-plusp (%bignum-0-or-plusp y (%bignum-length y)))
1798              (x (if x-plusp x (negate-bignum x nil)))
1799              (y (if y-plusp y (negate-bignum y nil)))
1800              (len-x (%bignum-length x))
1801              (len-y (%bignum-length y)))
1802         (multiple-value-bind (q r)
1803             (cond ((< len-y 2)
1804                    (bignum-truncate-single-digit x len-x y))
1805                   ((plusp (bignum-compare y x))
1806                    (let ((res (%allocate-bignum len-x)))
1807                      (dotimes (i len-x)
1808                        (setf (%bignum-ref res i) (%bignum-ref x i)))
1809                      (values 0 res)))
1810                   (t
1811                    (let ((len-x+1 (1+ len-x)))
1812                      (setf truncate-x (%allocate-bignum len-x+1))
1813                      (setf truncate-y (%allocate-bignum (1+ len-y)))
1814                      (let ((y-shift (shift-y-for-truncate y)))
1815                        (shift-and-store-truncate-buffers x len-x y
1816                                                          len-y y-shift)
1817                        (values (return-quotient-leaving-remainder len-x+1
1818                                                                   len-y)
1819                                ;; Now that RETURN-QUOTIENT-LEAVING-REMAINDER
1820                                ;; has executed, we just tidy up the remainder
1821                                ;; (in TRUNCATE-X) and return it.
1822                                (cond
1823                                  ((zerop y-shift)
1824                                   (let ((res (%allocate-bignum len-y)))
1825                                     (declare (type bignum-type res))
1826                                     (bignum-replace res truncate-x :end2 len-y)
1827                                     (%normalize-bignum res len-y)))
1828                                  (t
1829                                   (shift-right-unaligned
1830                                    truncate-x 0 y-shift len-y
1831                                    ((= j res-len-1)
1832                                     (setf (%bignum-ref res j)
1833                                           (%ashr (%bignum-ref truncate-x i)
1834                                                  y-shift))
1835                                     (%normalize-bignum res res-len))
1836                                    res))))))))
1837           (let ((quotient (cond ((eq x-plusp y-plusp) q)
1838                                 ((typep q 'fixnum) (the fixnum (- q)))
1839                                 (t (negate-bignum-in-place q))))
1840                 (rem (cond (x-plusp r)
1841                            ((typep r 'fixnum) (the fixnum (- r)))
1842                            (t (negate-bignum-in-place r)))))
1843             (values (if (typep quotient 'fixnum)
1844                         quotient
1845                         (%normalize-bignum quotient (%bignum-length quotient)))
1846                     (if (typep rem 'fixnum)
1847                         rem
1848                         (%normalize-bignum rem (%bignum-length rem))))))))))
1849
1850 \f
1851 ;;;; There used to be a pile of code for implementing division for bignum digits
1852 ;;;; for machines that don't have a 2*digit-size by digit-size divide instruction.
1853 ;;;; This happens to be most machines, but all the SBCL ports seem to be content
1854 ;;;; to implement SB-BIGNUM:%FLOOR as a VOP rather than using the code here.
1855 ;;;; So it's been deleted.  --njf, 2007-02-04
1856 \f
1857 ;;;; general utilities
1858
1859 ;;; Internal in-place operations use this to fixup remaining digits in the
1860 ;;; incoming data, such as in-place shifting. This is basically the same as
1861 ;;; the first form in %NORMALIZE-BIGNUM, but we return the length of the buffer
1862 ;;; instead of shrinking the bignum.
1863 #!-sb-fluid (declaim (sb!ext:maybe-inline %normalize-bignum-buffer))
1864 (defun %normalize-bignum-buffer (result len)
1865   (declare (type bignum-type result)
1866            (type bignum-index len))
1867   (unless (= len 1)
1868     (do ((next-digit (%bignum-ref result (- len 2))
1869                      (%bignum-ref result (- len 2)))
1870          (sign-digit (%bignum-ref result (1- len)) next-digit))
1871         ((not (zerop (logxor sign-digit (%ashr next-digit (1- digit-size))))))
1872         (decf len)
1873         (setf (%bignum-ref result len) 0)
1874         (when (= len 1)
1875               (return))))
1876   len)
1877
1878 ;;; This drops the last digit if it is unnecessary sign information. It repeats
1879 ;;; this as needed, possibly ending with a fixnum. If the resulting length from
1880 ;;; shrinking is one, see whether our one word is a fixnum. Shift the possible
1881 ;;; fixnum bits completely out of the word, and compare this with shifting the
1882 ;;; sign bit all the way through. If the bits are all 1's or 0's in both words,
1883 ;;; then there are just sign bits between the fixnum bits and the sign bit. If
1884 ;;; we do have a fixnum, shift it over for the two low-tag bits.
1885 (defun %normalize-bignum (result len)
1886   (declare (type bignum-type result)
1887            (type bignum-index len)
1888            (inline %normalize-bignum-buffer))
1889   (let ((newlen (%normalize-bignum-buffer result len)))
1890     (declare (type bignum-index newlen))
1891     (unless (= newlen len)
1892       (%bignum-set-length result newlen))
1893     (if (= newlen 1)
1894         (let ((digit (%bignum-ref result 0)))
1895           (if (= (%ashr digit sb!vm:n-positive-fixnum-bits)
1896                  (%ashr digit (1- digit-size)))
1897               (%fixnum-digit-with-correct-sign digit)
1898               result))
1899         result)))
1900
1901 ;;; This drops the last digit if it is unnecessary sign information. It
1902 ;;; repeats this as needed, possibly ending with a fixnum magnitude but never
1903 ;;; returning a fixnum.
1904 (defun %mostly-normalize-bignum (result len)
1905   (declare (type bignum-type result)
1906            (type bignum-index len)
1907            (inline %normalize-bignum-buffer))
1908   (let ((newlen (%normalize-bignum-buffer result len)))
1909     (declare (type bignum-index newlen))
1910     (unless (= newlen len)
1911       (%bignum-set-length result newlen))
1912     result))
1913 \f
1914 ;;;; hashing
1915
1916 ;;; the bignum case of the SXHASH function
1917 (defun sxhash-bignum (x)
1918   (let ((result 316495330))
1919     (declare (type fixnum result))
1920     (dotimes (i (%bignum-length x))
1921       (declare (type index i))
1922       (let ((xi (%bignum-ref x i)))
1923         (mixf result
1924               (logand most-positive-fixnum
1925                       xi
1926                       (ash xi -7)))))
1927     result))