src/code/numbers.lisp

   1 ;;;; This file contains the definitions of most number functions.
   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!KERNEL")
  13 \f
  14 ;;;; the NUMBER-DISPATCH macro
  15
  16 (eval-when (:compile-toplevel :load-toplevel :execute)
  17
  18 ;;; Grovel an individual case to NUMBER-DISPATCH, augmenting RESULT
  19 ;;; with the type dispatches and bodies. Result is a tree built of
  20 ;;; alists representing the dispatching off each arg (in order). The
  21 ;;; leaf is the body to be executed in that case.
  22 (defun parse-number-dispatch (vars result types var-types body)
  23   (cond ((null vars)
  24          (unless (null types) (error "More types than vars."))
  25          (when (cdr result)
  26            (error "Duplicate case: ~S." body))
  27          (setf (cdr result)
  28                (sublis var-types body :test #'equal)))
  29         ((null types)
  30          (error "More vars than types."))
  31         (t
  32          (flet ((frob (var type)
  33                   (parse-number-dispatch
  34                    (rest vars)
  35                    (or (assoc type (cdr result) :test #'equal)
  36                        (car (setf (cdr result)
  37                                   (acons type nil (cdr result)))))
  38                    (rest types)
  39                    (acons `(dispatch-type ,var) type var-types)
  40                    body)))
  41            (let ((type (first types))
  42                  (var (first vars)))
  43              (if (and (consp type) (eq (first type) 'foreach))
  44                  (dolist (type (rest type))
  45                    (frob var type))
  46                  (frob var type)))))))
  47
  48 ;;; our guess for the preferred order in which to do type tests
  49 ;;; (cheaper and/or more probable first.)
  50 (defparameter *type-test-ordering*
  51   '(fixnum single-float double-float integer #!+long-float long-float bignum
  52     complex ratio))
  53
  54 ;;; Should TYPE1 be tested before TYPE2?
  55 (defun type-test-order (type1 type2)
  56   (let ((o1 (position type1 *type-test-ordering*))
  57         (o2 (position type2 *type-test-ordering*)))
  58     (cond ((not o1) nil)
  59           ((not o2) t)
  60           (t
  61            (< o1 o2)))))
  62
  63 ;;; Return an ETYPECASE form that does the type dispatch, ordering the
  64 ;;; cases for efficiency.
  65 ;;; Check for some simple to detect problematic cases where the caller
  66 ;;; used types that are not disjoint and where this may lead to
  67 ;;; unexpected behaviour of the generated form, for example making
  68 ;;; a clause unreachable, and throw an error if such a case is found.
  69 ;;; An example:
  70 ;;;   (number-dispatch ((var1 integer) (var2 float))
  71 ;;;     ((fixnum single-float) a)
  72 ;;;     ((integer float) b))
  73 ;;; Even though the types are not reordered here, the generated form,
  74 ;;; basically
  75 ;;;   (etypecase var1
  76 ;;;     (fixnum (etypecase var2
  77 ;;;               (single-float a)))
  78 ;;;     (integer (etypecase var2
  79 ;;;                (float b))))
  80 ;;; would fail at runtime if given var1 fixnum and var2 double-float,
  81 ;;; even though the second clause matches this signature. To catch
  82 ;;; this earlier than runtime we throw an error already here.
  83 (defun generate-number-dispatch (vars error-tags cases)
  84   (if vars
  85       (let ((var (first vars))
  86             (cases (sort cases #'type-test-order :key #'car)))
  87         (flet ((error-if-sub-or-supertype (type1 type2)
  88                  (when (or (subtypep type1 type2)
  89                            (subtypep type2 type1))
  90                    (error "Types not disjoint: ~S ~S." type1 type2)))
  91                (error-if-supertype (type1 type2)
  92                  (when (subtypep type2 type1)
  93                    (error "Type ~S ordered before subtype ~S."
  94                           type1 type2)))
  95                (test-type-pairs (fun)
  96                  ;; Apply FUN to all (ordered) pairs of types from the
  97                  ;; cases.
  98                  (mapl (lambda (cases)
  99                          (when (cdr cases)
 100                            (let ((type1 (caar cases)))
 101                              (dolist (case (cdr cases))
 102                                (funcall fun type1 (car case))))))
 103                        cases)))
 104           ;; For the last variable throw an error if a type is followed
 105           ;; by a subtype, for all other variables additionally if a
 106           ;; type is followed by a supertype.
 107           (test-type-pairs (if (cdr vars)
 108                                #'error-if-sub-or-supertype
 109                                #'error-if-supertype)))
 110         `((typecase ,var
 111             ,@(mapcar (lambda (case)
 112                         `(,(first case)
 113                           ,@(generate-number-dispatch (rest vars)
 114                                                       (rest error-tags)
 115                                                       (cdr case))))
 116                       cases)
 117             (t (go ,(first error-tags))))))
 118       cases))
 119
 120 ) ; EVAL-WHEN
 121
 122 ;;; This is a vaguely case-like macro that does number cross-product
 123 ;;; dispatches. The Vars are the variables we are dispatching off of.
 124 ;;; The Type paired with each Var is used in the error message when no
 125 ;;; case matches. Each case specifies a Type for each var, and is
 126 ;;; executed when that signature holds. A type may be a list
 127 ;;; (FOREACH Each-Type*), causing that case to be repeatedly
 128 ;;; instantiated for every Each-Type. In the body of each case, any
 129 ;;; list of the form (DISPATCH-TYPE Var-Name) is substituted with the
 130 ;;; type of that var in that instance of the case.
 131 ;;;
 132 ;;; As an alternate to a case spec, there may be a form whose CAR is a
 133 ;;; symbol. In this case, we apply the CAR of the form to the CDR and
 134 ;;; treat the result of the call as a list of cases. This process is
 135 ;;; not applied recursively.
 136 ;;;
 137 ;;; Be careful when using non-disjoint types in different cases for the
 138 ;;; same variable. Some uses will behave as intended, others not, as the
 139 ;;; variables are dispatched off sequentially and clauses are reordered
 140 ;;; for efficiency. Some, but not all, problematic cases are detected
 141 ;;; and lead to a compile time error; see GENERATE-NUMBER-DISPATCH above
 142 ;;; for an example.
 143 (defmacro number-dispatch (var-specs &body cases)
 144   (let ((res (list nil))
 145         (vars (mapcar #'car var-specs))
 146         (block (gensym)))
 147     (dolist (case cases)
 148       (if (symbolp (first case))
 149           (let ((cases (apply (symbol-function (first case)) (rest case))))
 150             (dolist (case cases)
 151               (parse-number-dispatch vars res (first case) nil (rest case))))
 152           (parse-number-dispatch vars res (first case) nil (rest case))))
 153
 154     (collect ((errors)
 155               (error-tags))
 156       (dolist (spec var-specs)
 157         (let ((var (first spec))
 158               (type (second spec))
 159               (tag (gensym)))
 160           (error-tags tag)
 161           (errors tag)
 162           (errors `(return-from
 163                     ,block
 164                     (error 'simple-type-error :datum ,var
 165                            :expected-type ',type
 166                            :format-control
 167                            "~@<Argument ~A is not a ~S: ~2I~_~S~:>"
 168                            :format-arguments
 169                            (list ',var ',type ,var))))))
 170
 171       `(block ,block
 172          (tagbody
 173            (return-from ,block
 174                         ,@(generate-number-dispatch vars (error-tags)
 175                                                     (cdr res)))
 176            ,@(errors))))))
 177 \f
 178 ;;;; binary operation dispatching utilities
 179
 180 (eval-when (:compile-toplevel :execute)
 181
 182 ;;; Return NUMBER-DISPATCH forms for rational X float.
 183 (defun float-contagion (op x y &optional (rat-types '(fixnum bignum ratio)))
 184   `(((single-float single-float) (,op ,x ,y))
 185     (((foreach ,@rat-types)
 186       (foreach single-float double-float #!+long-float long-float))
 187      (,op (coerce ,x '(dispatch-type ,y)) ,y))
 188     (((foreach single-float double-float #!+long-float long-float)
 189       (foreach ,@rat-types))
 190      (,op ,x (coerce ,y '(dispatch-type ,x))))
 191     #!+long-float
 192     (((foreach single-float double-float long-float) long-float)
 193      (,op (coerce ,x 'long-float) ,y))
 194     #!+long-float
 195     ((long-float (foreach single-float double-float))
 196      (,op ,x (coerce ,y 'long-float)))
 197     (((foreach single-float double-float) double-float)
 198      (,op (coerce ,x 'double-float) ,y))
 199     ((double-float single-float)
 200      (,op ,x (coerce ,y 'double-float)))))
 201
 202 ;;; Return NUMBER-DISPATCH forms for bignum X fixnum.
 203 (defun bignum-cross-fixnum (fix-op big-op)
 204   `(((fixnum fixnum) (,fix-op x y))
 205     ((fixnum bignum)
 206      (,big-op (make-small-bignum x) y))
 207     ((bignum fixnum)
 208      (,big-op x (make-small-bignum y)))
 209     ((bignum bignum)
 210      (,big-op x y))))
 211
 212 ) ; EVAL-WHEN
 213 \f
 214 ;;;; canonicalization utilities
 215
 216 ;;; If IMAGPART is 0, return REALPART, otherwise make a complex. This is
 217 ;;; used when we know that REALPART and IMAGPART are the same type, but
 218 ;;; rational canonicalization might still need to be done.
 219 #!-sb-fluid (declaim (inline canonical-complex))
 220 (defun canonical-complex (realpart imagpart)
 221   (if (eql imagpart 0)
 222       realpart
 223       (cond #!+long-float
 224             ((and (typep realpart 'long-float)
 225                   (typep imagpart 'long-float))
 226              (truly-the (complex long-float) (complex realpart imagpart)))
 227             ((and (typep realpart 'double-float)
 228                   (typep imagpart 'double-float))
 229              (truly-the (complex double-float) (complex realpart imagpart)))
 230             ((and (typep realpart 'single-float)
 231                   (typep imagpart 'single-float))
 232              (truly-the (complex single-float) (complex realpart imagpart)))
 233             (t
 234              (%make-complex realpart imagpart)))))
 235
 236 ;;; Given a numerator and denominator with the GCD already divided
 237 ;;; out, make a canonical rational. We make the denominator positive,
 238 ;;; and check whether it is 1.
 239 #!-sb-fluid (declaim (inline build-ratio))
 240 (defun build-ratio (num den)
 241   (multiple-value-bind (num den)
 242       (if (minusp den)
 243           (values (- num) (- den))
 244           (values num den))
 245     (cond
 246       ((eql den 0)
 247        (error 'division-by-zero
 248               :operands (list num den)
 249               :operation 'build-ratio))
 250       ((eql den 1) num)
 251       (t (%make-ratio num den)))))
 252
 253 ;;; Truncate X and Y, but bum the case where Y is 1.
 254 #!-sb-fluid (declaim (inline maybe-truncate))
 255 (defun maybe-truncate (x y)
 256   (if (eql y 1)
 257       x
 258       (truncate x y)))
 259 \f
 260 ;;;; COMPLEXes
 261
 262 (defun complex (realpart &optional (imagpart 0))
 263   #!+sb-doc
 264   "Return a complex number with the specified real and imaginary components."
 265   (declare (explicit-check))
 266   (flet ((%%make-complex (realpart imagpart)
 267            (cond #!+long-float
 268                  ((and (typep realpart 'long-float)
 269                        (typep imagpart 'long-float))
 270                   (truly-the (complex long-float)
 271                              (complex realpart imagpart)))
 272                  ((and (typep realpart 'double-float)
 273                        (typep imagpart 'double-float))
 274                   (truly-the (complex double-float)
 275                              (complex realpart imagpart)))
 276                  ((and (typep realpart 'single-float)
 277                        (typep imagpart 'single-float))
 278                   (truly-the (complex single-float)
 279                              (complex realpart imagpart)))
 280                  (t
 281                   (%make-complex realpart imagpart)))))
 282   (number-dispatch ((realpart real) (imagpart real))
 283     ((rational rational)
 284      (canonical-complex realpart imagpart))
 285     (float-contagion %%make-complex realpart imagpart (rational)))))
 286
 287 (defun realpart (number)
 288   #!+sb-doc
 289   "Extract the real part of a number."
 290   (etypecase number
 291     #!+long-float
 292     ((complex long-float)
 293      (truly-the long-float (realpart number)))
 294     ((complex double-float)
 295      (truly-the double-float (realpart number)))
 296     ((complex single-float)
 297      (truly-the single-float (realpart number)))
 298     ((complex rational)
 299      (%realpart number))
 300     (number
 301      number)))
 302
 303 (defun imagpart (number)
 304   #!+sb-doc
 305   "Extract the imaginary part of a number."
 306   (etypecase number
 307     #!+long-float
 308     ((complex long-float)
 309      (truly-the long-float (imagpart number)))
 310     ((complex double-float)
 311      (truly-the double-float (imagpart number)))
 312     ((complex single-float)
 313      (truly-the single-float (imagpart number)))
 314     ((complex rational)
 315      (%imagpart number))
 316     (float
 317      (* 0 number))
 318     (number
 319      0)))
 320
 321 (defun conjugate (number)
 322   #!+sb-doc
 323   "Return the complex conjugate of NUMBER. For non-complex numbers, this is
 324   an identity."
 325   (declare (type number number) (explicit-check))
 326   (if (complexp number)
 327       (complex (realpart number) (- (imagpart number)))
 328       number))
 329
 330 (defun signum (number)
 331   #!+sb-doc
 332   "If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))."
 333   (declare (explicit-check))
 334   (if (zerop number)
 335       number
 336       (if (rationalp number)
 337           (if (plusp number) 1 -1)
 338           (/ number (abs number)))))
 339 \f
 340 ;;;; ratios
 341
 342 (defun numerator (number)
 343   #!+sb-doc
 344   "Return the numerator of NUMBER, which must be rational."
 345   (numerator number))
 346
 347 (defun denominator (number)
 348   #!+sb-doc
 349   "Return the denominator of NUMBER, which must be rational."
 350   (denominator number))
 351 \f
 352 ;;;; arithmetic operations
 353 ;;;;
 354 ;;;; IMPORTANT NOTE: Accessing &REST arguments with NTH is actually extremely
 355 ;;;; efficient in SBCL, as is taking their LENGTH -- so this code is very
 356 ;;;; clever instead of being charmingly naive. Please check that "obvious"
 357 ;;;; improvements don't actually ruin performance.
 358 ;;;;
 359 ;;;; (Granted that the difference between very clever and charmingly naivve
 360 ;;;; can sometimes be sliced exceedingly thing...)
 361
 362 (macrolet ((define-arith (op init doc)
 363              #!-sb-doc (declare (ignore doc))
 364              `(defun ,op (&rest numbers)
 365                 (declare (explicit-check))
 366                 #!+sb-doc ,doc
 367                 (if numbers
 368                     (let ((result (the number (fast-&rest-nth 0 numbers))))
 369                       (do-rest-arg ((n) numbers 1 result)
 370                         (setq result (,op result n))))
 371                     ,init))))
 372   (define-arith + 0
 373     "Return the sum of its arguments. With no args, returns 0.")
 374   (define-arith * 1
 375     "Return the product of its arguments. With no args, returns 1."))
 376
 377 (defun - (number &rest more-numbers)
 378   #!+sb-doc
 379   "Subtract the second and all subsequent arguments from the first;
 380   or with one argument, negate the first argument."
 381   (declare (explicit-check))
 382   (if more-numbers
 383       (let ((result number))
 384         (do-rest-arg ((n) more-numbers 0 result)
 385           (setf result (- result n))))
 386       (- number)))
 387
 388 (defun / (number &rest more-numbers)
 389   #!+sb-doc
 390   "Divide the first argument by each of the following arguments, in turn.
 391   With one argument, return reciprocal."
 392   (declare (explicit-check))
 393   (if more-numbers
 394       (let ((result number))
 395         (do-rest-arg ((n) more-numbers 0 result)
 396           (setf result (/ result n))))
 397       (/ number)))
 398
 399 (defun 1+ (number)
 400   #!+sb-doc
 401   "Return NUMBER + 1."
 402   (declare (explicit-check))
 403   (1+ number))
 404
 405 (defun 1- (number)
 406   #!+sb-doc
 407   "Return NUMBER - 1."
 408   (declare (explicit-check))
 409   (1- number))
 410
 411 (eval-when (:compile-toplevel)
 412
 413 (sb!xc:defmacro two-arg-+/- (name op big-op)
 414   `(defun ,name (x y)
 415      (number-dispatch ((x number) (y number))
 416        (bignum-cross-fixnum ,op ,big-op)
 417        (float-contagion ,op x y)
 418
 419        ((complex complex)
 420         (canonical-complex (,op (realpart x) (realpart y))
 421                            (,op (imagpart x) (imagpart y))))
 422        (((foreach bignum fixnum ratio single-float double-float
 423                   #!+long-float long-float) complex)
 424         (complex (,op x (realpart y)) (,op 0 (imagpart y))))
 425        ((complex (or rational float))
 426         (complex (,op (realpart x) y) (,op (imagpart x) 0)))
 427
 428        (((foreach fixnum bignum) ratio)
 429         (let* ((dy (denominator y))
 430                (n (,op (* x dy) (numerator y))))
 431           (%make-ratio n dy)))
 432        ((ratio integer)
 433         (let* ((dx (denominator x))
 434                (n (,op (numerator x) (* y dx))))
 435           (%make-ratio n dx)))
 436        ((ratio ratio)
 437         (let* ((nx (numerator x))
 438                (dx (denominator x))
 439                (ny (numerator y))
 440                (dy (denominator y))
 441                (g1 (gcd dx dy)))
 442           (if (eql g1 1)
 443               (%make-ratio (,op (* nx dy) (* dx ny)) (* dx dy))
 444               (let* ((t1 (,op (* nx (truncate dy g1)) (* (truncate dx g1) ny)))
 445                      (g2 (gcd t1 g1))
 446                      (t2 (truncate dx g1)))
 447                 (cond ((eql t1 0) 0)
 448                       ((eql g2 1)
 449                        (%make-ratio t1 (* t2 dy)))
 450                       (t (let* ((nn (truncate t1 g2))
 451                                 (t3 (truncate dy g2))
 452                                 (nd (if (eql t2 1) t3 (* t2 t3))))
 453                            (if (eql nd 1) nn (%make-ratio nn nd))))))))))))
 454
 455 ) ; EVAL-WHEN
 456
 457 (two-arg-+/- two-arg-+ + add-bignums)
 458 (two-arg-+/- two-arg-- - subtract-bignum)
 459
 460 (defun two-arg-* (x y)
 461   (flet ((integer*ratio (x y)
 462            (if (eql x 0) 0
 463                (let* ((ny (numerator y))
 464                       (dy (denominator y))
 465                       (gcd (gcd x dy)))
 466                  (if (eql gcd 1)
 467                      (%make-ratio (* x ny) dy)
 468                      (let ((nn (* (truncate x gcd) ny))
 469                            (nd (truncate dy gcd)))
 470                        (if (eql nd 1)
 471                            nn
 472                            (%make-ratio nn nd)))))))
 473          (complex*real (x y)
 474            (canonical-complex (* (realpart x) y) (* (imagpart x) y))))
 475     (number-dispatch ((x number) (y number))
 476       (float-contagion * x y)
 477
 478       ((fixnum fixnum) (multiply-fixnums x y))
 479       ((bignum fixnum) (multiply-bignum-and-fixnum x y))
 480       ((fixnum bignum) (multiply-bignum-and-fixnum y x))
 481       ((bignum bignum) (multiply-bignums x y))
 482
 483       ((complex complex)
 484        (let* ((rx (realpart x))
 485               (ix (imagpart x))
 486               (ry (realpart y))
 487               (iy (imagpart y)))
 488          (canonical-complex (- (* rx ry) (* ix iy)) (+ (* rx iy) (* ix ry)))))
 489       (((foreach bignum fixnum ratio single-float double-float
 490                  #!+long-float long-float)
 491         complex)
 492        (complex*real y x))
 493       ((complex (or rational float))
 494        (complex*real x y))
 495
 496       (((foreach bignum fixnum) ratio) (integer*ratio x y))
 497       ((ratio integer) (integer*ratio y x))
 498       ((ratio ratio)
 499        (let* ((nx (numerator x))
 500               (dx (denominator x))
 501               (ny (numerator y))
 502               (dy (denominator y))
 503               (g1 (gcd nx dy))
 504               (g2 (gcd dx ny)))
 505          (build-ratio (* (maybe-truncate nx g1)
 506                          (maybe-truncate ny g2))
 507                       (* (maybe-truncate dx g2)
 508                          (maybe-truncate dy g1))))))))
 509
 510 ;;; Divide two integers, producing a canonical rational. If a fixnum,
 511 ;;; we see whether they divide evenly before trying the GCD. In the
 512 ;;; bignum case, we don't bother, since bignum division is expensive,
 513 ;;; and the test is not very likely to succeed.
 514 (defun integer-/-integer (x y)
 515   (if (and (typep x 'fixnum) (typep y 'fixnum))
 516       (multiple-value-bind (quo rem) (truncate x y)
 517         (if (zerop rem)
 518             quo
 519             (let ((gcd (gcd x y)))
 520               (declare (fixnum gcd))
 521               (if (eql gcd 1)
 522                   (build-ratio x y)
 523                   (build-ratio (truncate x gcd) (truncate y gcd))))))
 524       (let ((gcd (gcd x y)))
 525         (if (eql gcd 1)
 526             (build-ratio x y)
 527             (build-ratio (truncate x gcd) (truncate y gcd))))))
 528
 529 (defun two-arg-/ (x y)
 530   (number-dispatch ((x number) (y number))
 531     (float-contagion / x y (ratio integer))
 532
 533     ((complex complex)
 534      (let* ((rx (realpart x))
 535             (ix (imagpart x))
 536             (ry (realpart y))
 537             (iy (imagpart y)))
 538        (if (> (abs ry) (abs iy))
 539            (let* ((r (/ iy ry))
 540                   (dn (* ry (+ 1 (* r r)))))
 541              (canonical-complex (/ (+ rx (* ix r)) dn)
 542                                 (/ (- ix (* rx r)) dn)))
 543            (let* ((r (/ ry iy))
 544                   (dn (* iy (+ 1 (* r r)))))
 545              (canonical-complex (/ (+ (* rx r) ix) dn)
 546                                 (/ (- (* ix r) rx) dn))))))
 547     (((foreach integer ratio single-float double-float) complex)
 548      (let* ((ry (realpart y))
 549             (iy (imagpart y)))
 550        (if (> (abs ry) (abs iy))
 551            (let* ((r (/ iy ry))
 552                   (dn (* ry (+ 1 (* r r)))))
 553              (canonical-complex (/ x dn)
 554                                 (/ (- (* x r)) dn)))
 555            (let* ((r (/ ry iy))
 556                   (dn (* iy (+ 1 (* r r)))))
 557              (canonical-complex (/ (* x r) dn)
 558                                 (/ (- x) dn))))))
 559     ((complex (or rational float))
 560      (canonical-complex (/ (realpart x) y)
 561                         (/ (imagpart x) y)))
 562
 563     ((ratio ratio)
 564      (let* ((nx (numerator x))
 565             (dx (denominator x))
 566             (ny (numerator y))
 567             (dy (denominator y))
 568             (g1 (gcd nx ny))
 569             (g2 (gcd dx dy)))
 570        (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2))
 571                     (* (maybe-truncate dx g2) (maybe-truncate ny g1)))))
 572
 573     ((integer integer)
 574      (integer-/-integer x y))
 575
 576     ((integer ratio)
 577      (if (zerop x)
 578          0
 579          (let* ((ny (numerator y))
 580                 (dy (denominator y))
 581                 (gcd (gcd x ny)))
 582            (build-ratio (* (maybe-truncate x gcd) dy)
 583                         (maybe-truncate ny gcd)))))
 584
 585     ((ratio integer)
 586      (let* ((nx (numerator x))
 587             (gcd (gcd nx y)))
 588        (build-ratio (maybe-truncate nx gcd)
 589                     (* (maybe-truncate y gcd) (denominator x)))))))
 590
 591 (defun %negate (n)
 592   (declare (explicit-check))
 593   (number-dispatch ((n number))
 594     (((foreach fixnum single-float double-float #!+long-float long-float))
 595      (%negate n))
 596     ((bignum)
 597      (negate-bignum n))
 598     ((ratio)
 599      (%make-ratio (- (numerator n)) (denominator n)))
 600     ((complex)
 601      (complex (- (realpart n)) (- (imagpart n))))))
 602 \f
 603 ;;;; TRUNCATE and friends
 604
 605 (defun truncate (number &optional (divisor 1))
 606   #!+sb-doc
 607   "Return number (or number/divisor) as an integer, rounded toward 0.
 608   The second returned value is the remainder."
 609   (declare (explicit-check))
 610   (macrolet ((truncate-float (rtype)
 611                `(let* ((float-div (coerce divisor ',rtype))
 612                        (res (%unary-truncate (/ number float-div))))
 613                   (values res
 614                           (- number
 615                              (* (coerce res ',rtype) float-div))))))
 616     (number-dispatch ((number real) (divisor real))
 617       ((fixnum fixnum) (truncate number divisor))
 618       (((foreach fixnum bignum) ratio)
 619        (if (= (numerator divisor) 1)
 620            (values (* number (denominator divisor)) 0)
 621            (multiple-value-bind (quot rem)
 622                (truncate (* number (denominator divisor))
 623                          (numerator divisor))
 624              (values quot (/ rem (denominator divisor))))))
 625       ((fixnum bignum)
 626        (bignum-truncate (make-small-bignum number) divisor))
 627       ((ratio (or float rational))
 628        (let ((q (truncate (numerator number)
 629                           (* (denominator number) divisor))))
 630          (values q (- number (* q divisor)))))
 631       ((bignum fixnum)
 632        (bignum-truncate number (make-small-bignum divisor)))
 633       ((bignum bignum)
 634        (bignum-truncate number divisor))
 635
 636       (((foreach single-float double-float #!+long-float long-float)
 637         (or rational single-float))
 638        (if (eql divisor 1)
 639            (let ((res (%unary-truncate number)))
 640              (values res (- number (coerce res '(dispatch-type number)))))
 641            (truncate-float (dispatch-type number))))
 642       #!+long-float
 643       ((long-float (or single-float double-float long-float))
 644        (truncate-float long-float))
 645       #!+long-float
 646       (((foreach double-float single-float) long-float)
 647        (truncate-float long-float))
 648       ((double-float (or single-float double-float))
 649        (truncate-float double-float))
 650       ((single-float double-float)
 651        (truncate-float double-float))
 652       (((foreach fixnum bignum ratio)
 653         (foreach single-float double-float #!+long-float long-float))
 654        (truncate-float (dispatch-type divisor))))))
 655
 656 (defun %multiply-high (x y)
 657   (declare (type word x y))
 658   (%multiply-high x y))
 659
 660 (defun floor (number &optional (divisor 1))
 661   #!+sb-doc
 662   "Return the greatest integer not greater than number, or number/divisor.
 663   The second returned value is (mod number divisor)."
 664   (declare (explicit-check))
 665   (floor number divisor))
 666
 667 (defun ceiling (number &optional (divisor 1))
 668   #!+sb-doc
 669   "Return the smallest integer not less than number, or number/divisor.
 670   The second returned value is the remainder."
 671   (declare (explicit-check))
 672   (ceiling number divisor))
 673
 674 (defun rem (number divisor)
 675   #!+sb-doc
 676   "Return second result of TRUNCATE."
 677   (declare (explicit-check))
 678   (rem number divisor))
 679
 680 (defun mod (number divisor)
 681   #!+sb-doc
 682   "Return second result of FLOOR."
 683   (declare (explicit-check))
 684   (mod number divisor))
 685
 686 (defun round (number &optional (divisor 1))
 687   #!+sb-doc
 688   "Rounds number (or number/divisor) to nearest integer.
 689   The second returned value is the remainder."
 690   (declare (explicit-check))
 691   (if (eql divisor 1)
 692       (round number)
 693       (multiple-value-bind (tru rem) (truncate number divisor)
 694         (if (zerop rem)
 695             (values tru rem)
 696             (let ((thresh (/ (abs divisor) 2)))
 697               (cond ((or (> rem thresh)
 698                          (and (= rem thresh) (oddp tru)))
 699                      (if (minusp divisor)
 700                          (values (- tru 1) (+ rem divisor))
 701                          (values (+ tru 1) (- rem divisor))))
 702                     ((let ((-thresh (- thresh)))
 703                        (or (< rem -thresh)
 704                            (and (= rem -thresh) (oddp tru))))
 705                      (if (minusp divisor)
 706                          (values (+ tru 1) (- rem divisor))
 707                          (values (- tru 1) (+ rem divisor))))
 708                     (t (values tru rem))))))))
 709
 710 (defmacro !define-float-rounding-function (name op doc)
 711   `(defun ,name (number &optional (divisor 1))
 712     ,doc
 713     (multiple-value-bind (res rem) (,op number divisor)
 714       (values (float res (if (floatp rem) rem 1.0)) rem))))
 715
 716 ;;; Declare these guys inline to let them get optimized a little.
 717 ;;; ROUND and FROUND are not declared inline since they seem too
 718 ;;; obscure and too big to inline-expand by default. Also, this gives
 719 ;;; the compiler a chance to pick off the unary float case.
 720 #!-sb-fluid (declaim (inline fceiling ffloor ftruncate))
 721 (defun ftruncate (number &optional (divisor 1))
 722   #!+sb-doc
 723   "Same as TRUNCATE, but returns first value as a float."
 724   (declare (explicit-check))
 725   (macrolet ((ftruncate-float (rtype)
 726                `(let* ((float-div (coerce divisor ',rtype))
 727                        (res (%unary-ftruncate (/ number float-div))))
 728                   (values res
 729                           (- number
 730                              (* (coerce res ',rtype) float-div))))))
 731     (number-dispatch ((number real) (divisor real))
 732       (((foreach fixnum bignum ratio) (or fixnum bignum ratio))
 733        (multiple-value-bind (q r)
 734            (truncate number divisor)
 735          (values (float q) r)))
 736       (((foreach single-float double-float #!+long-float long-float)
 737         (or rational single-float))
 738        (if (eql divisor 1)
 739            (let ((res (%unary-ftruncate number)))
 740              (values res (- number (coerce res '(dispatch-type number)))))
 741            (ftruncate-float (dispatch-type number))))
 742       #!+long-float
 743       ((long-float (or single-float double-float long-float))
 744        (ftruncate-float long-float))
 745       #!+long-float
 746       (((foreach double-float single-float) long-float)
 747        (ftruncate-float long-float))
 748       ((double-float (or single-float double-float))
 749        (ftruncate-float double-float))
 750       ((single-float double-float)
 751        (ftruncate-float double-float))
 752       (((foreach fixnum bignum ratio)
 753         (foreach single-float double-float #!+long-float long-float))
 754        (ftruncate-float (dispatch-type divisor))))))
 755
 756 (defun ffloor (number &optional (divisor 1))
 757   #!+sb-doc
 758   "Same as FLOOR, but returns first value as a float."
 759   (declare (explicit-check))
 760   (multiple-value-bind (tru rem) (ftruncate number divisor)
 761     (if (and (not (zerop rem))
 762              (if (minusp divisor)
 763                  (plusp number)
 764                  (minusp number)))
 765         (values (1- tru) (+ rem divisor))
 766         (values tru rem))))
 767
 768 (defun fceiling (number &optional (divisor 1))
 769   #!+sb-doc
 770   "Same as CEILING, but returns first value as a float."
 771   (declare (explicit-check))
 772   (multiple-value-bind (tru rem) (ftruncate number divisor)
 773     (if (and (not (zerop rem))
 774              (if (minusp divisor)
 775                  (minusp number)
 776                  (plusp number)))
 777         (values (+ tru 1) (- rem divisor))
 778         (values tru rem))))
 779
 780 ;;; FIXME: this probably needs treatment similar to the use of
 781 ;;; %UNARY-FTRUNCATE for FTRUNCATE.
 782 (defun fround (number &optional (divisor 1))
 783   #!+sb-doc
 784   "Same as ROUND, but returns first value as a float."
 785   (declare (explicit-check))
 786   (multiple-value-bind (res rem)
 787       (round number divisor)
 788     (values (float res (if (floatp rem) rem 1.0)) rem)))
 789 \f
 790 ;;;; comparisons
 791
 792 (defun = (number &rest more-numbers)
 793   #!+sb-doc
 794   "Return T if all of its arguments are numerically equal, NIL otherwise."
 795   (declare (number number) (explicit-check))
 796   (do-rest-arg ((n i) more-numbers 0 t)
 797     (unless (= number n)
 798       (return (do-rest-arg ((n) more-numbers (1+ i))
 799                 (the number n)))))) ; for effect
 800
 801 (defun /= (number &rest more-numbers)
 802   #!+sb-doc
 803   "Return T if no two of its arguments are numerically equal, NIL otherwise."
 804   (declare (number number) (explicit-check))
 805   (if more-numbers
 806       (do ((n number (nth i more-numbers))
 807             (i 0 (1+ i)))
 808           ((>= i (length more-numbers))
 809            t)
 810         (do-rest-arg ((n2) more-numbers i)
 811           (when (= n n2)
 812             (return-from /= nil))))
 813       t))
 814
 815 (macrolet ((def (op doc)
 816              (declare (ignorable doc))
 817              `(defun ,op (number &rest more-numbers)
 818                 #!+sb-doc ,doc
 819                 (declare (explicit-check))
 820                 (let ((n1 number))
 821                   (declare (real n1))
 822                   (do-rest-arg ((n2 i) more-numbers 0 t)
 823                     (if (,op n1 n2)
 824                         (setf n1 n2)
 825                         (return (do-rest-arg ((n) more-numbers (1+ i))
 826                                   (the real n))))))))) ; for effect
 827   (def <  "Return T if its arguments are in strictly increasing order, NIL otherwise.")
 828   (def >  "Return T if its arguments are in strictly decreasing order, NIL otherwise.")
 829   (def <= "Return T if arguments are in strictly non-decreasing order, NIL otherwise.")
 830   (def >= "Return T if arguments are in strictly non-increasing order, NIL otherwise."))
 831
 832 (defun max (number &rest more-numbers)
 833   #!+sb-doc
 834   "Return the greatest of its arguments; among EQUALP greatest, return
 835 the first."
 836   (declare (explicit-check))
 837   (let ((n number))
 838     (declare (real n))
 839     (do-rest-arg ((arg) more-numbers 0 n)
 840       (when (> arg n)
 841         (setf n arg)))))
 842
 843 (defun min (number &rest more-numbers)
 844   #!+sb-doc
 845   "Return the least of its arguments; among EQUALP least, return
 846 the first."
 847   (declare (explicit-check))
 848   (let ((n number))
 849     (declare (real n))
 850     (do-rest-arg ((arg) more-numbers 0 n)
 851       (when (< arg n)
 852         (setf n arg)))))
 853
 854 (eval-when (:compile-toplevel :execute)
 855
 856 ;;; The INFINITE-X-FINITE-Y and INFINITE-Y-FINITE-X args tell us how
 857 ;;; to handle the case when X or Y is a floating-point infinity and
 858 ;;; the other arg is a rational. (Section 12.1.4.1 of the ANSI spec
 859 ;;; says that comparisons are done by converting the float to a
 860 ;;; rational when comparing with a rational, but infinities can't be
 861 ;;; converted to a rational, so we show some initiative and do it this
 862 ;;; way instead.)
 863 (defun basic-compare (op &key infinite-x-finite-y infinite-y-finite-x)
 864   `(((fixnum fixnum) (,op x y))
 865
 866     ((single-float single-float) (,op x y))
 867     #!+long-float
 868     (((foreach single-float double-float long-float) long-float)
 869      (,op (coerce x 'long-float) y))
 870     #!+long-float
 871     ((long-float (foreach single-float double-float))
 872      (,op x (coerce y 'long-float)))
 873     ((fixnum (foreach single-float double-float))
 874      (if (float-infinity-p y)
 875          ,infinite-y-finite-x
 876          ;; If the fixnum has an exact float representation, do a
 877          ;; float comparison. Otherwise do the slow float -> ratio
 878          ;; conversion.
 879          (multiple-value-bind (lo hi)
 880              (case '(dispatch-type y)
 881                (single-float
 882                 (values most-negative-exactly-single-float-fixnum
 883                         most-positive-exactly-single-float-fixnum))
 884                (double-float
 885                 (values most-negative-exactly-double-float-fixnum
 886                         most-positive-exactly-double-float-fixnum)))
 887            (if (<= lo y hi)
 888                (,op (coerce x '(dispatch-type y)) y)
 889                (,op x (rational y))))))
 890     (((foreach single-float double-float) fixnum)
 891      (if (eql y 0)
 892          (,op x (coerce 0 '(dispatch-type x)))
 893          (if (float-infinity-p x)
 894              ,infinite-x-finite-y
 895              ;; Likewise
 896              (multiple-value-bind (lo hi)
 897                  (case '(dispatch-type x)
 898                    (single-float
 899                     (values most-negative-exactly-single-float-fixnum
 900                             most-positive-exactly-single-float-fixnum))
 901                    (double-float
 902                     (values most-negative-exactly-double-float-fixnum
 903                             most-positive-exactly-double-float-fixnum)))
 904                (if (<= lo y hi)
 905                    (,op x (coerce y '(dispatch-type x)))
 906                    (,op (rational x) y))))))
 907     (((foreach single-float double-float) double-float)
 908      (,op (coerce x 'double-float) y))
 909     ((double-float single-float)
 910      (,op x (coerce y 'double-float)))
 911     (((foreach single-float double-float #!+long-float long-float) rational)
 912      (if (eql y 0)
 913          (,op x (coerce 0 '(dispatch-type x)))
 914          (if (float-infinity-p x)
 915              ,infinite-x-finite-y
 916              (,op (rational x) y))))
 917     (((foreach bignum fixnum ratio) float)
 918      (if (float-infinity-p y)
 919          ,infinite-y-finite-x
 920          (,op x (rational y))))))
 921 ) ; EVAL-WHEN
 922
 923 (macrolet ((def-two-arg-</> (name op ratio-arg1 ratio-arg2 &rest cases)
 924              `(defun ,name (x y)
 925                 (number-dispatch ((x real) (y real))
 926                                  (basic-compare
 927                                   ,op
 928                                   :infinite-x-finite-y
 929                                   (,op x (coerce 0 '(dispatch-type x)))
 930                                   :infinite-y-finite-x
 931                                   (,op (coerce 0 '(dispatch-type y)) y))
 932                                  (((foreach fixnum bignum) ratio)
 933                                   (,op x (,ratio-arg2 (numerator y)
 934                                                       (denominator y))))
 935                                  ((ratio integer)
 936                                   (,op (,ratio-arg1 (numerator x)
 937                                                     (denominator x))
 938                                        y))
 939                                  ((ratio ratio)
 940                                   (,op (* (numerator   (truly-the ratio x))
 941                                           (denominator (truly-the ratio y)))
 942                                        (* (numerator   (truly-the ratio y))
 943                                           (denominator (truly-the ratio x)))))
 944                                  ,@cases))))
 945   (def-two-arg-</> two-arg-< < floor ceiling
 946     ((fixnum bignum)
 947      (bignum-plus-p y))
 948     ((bignum fixnum)
 949      (not (bignum-plus-p x)))
 950     ((bignum bignum)
 951      (minusp (bignum-compare x y))))
 952   (def-two-arg-</> two-arg-> > ceiling floor
 953     ((fixnum bignum)
 954      (not (bignum-plus-p y)))
 955     ((bignum fixnum)
 956      (bignum-plus-p x))
 957     ((bignum bignum)
 958      (plusp (bignum-compare x y)))))
 959
 960 (defun two-arg-= (x y)
 961   (number-dispatch ((x number) (y number))
 962     (basic-compare =
 963                    ;; An infinite value is never equal to a finite value.
 964                    :infinite-x-finite-y nil
 965                    :infinite-y-finite-x nil)
 966     ((fixnum (or bignum ratio)) nil)
 967
 968     ((bignum (or fixnum ratio)) nil)
 969     ((bignum bignum)
 970      (zerop (bignum-compare x y)))
 971
 972     ((ratio integer) nil)
 973     ((ratio ratio)
 974      (and (eql (numerator x) (numerator y))
 975           (eql (denominator x) (denominator y))))
 976
 977     ((complex complex)
 978      (and (= (realpart x) (realpart y))
 979           (= (imagpart x) (imagpart y))))
 980     (((foreach fixnum bignum ratio single-float double-float
 981                #!+long-float long-float) complex)
 982      (and (= x (realpart y))
 983           (zerop (imagpart y))))
 984     ((complex (or float rational))
 985      (and (= (realpart x) y)
 986           (zerop (imagpart x))))))
 987 \f
 988 ;;;; logicals
 989
 990 (macrolet ((def (op init doc)
 991              #!-sb-doc (declare (ignore doc))
 992              `(defun ,op (&rest integers)
 993                 #!+sb-doc ,doc
 994                 (declare (explicit-check))
 995                 (if integers
 996                     (do ((result (fast-&rest-nth 0 integers)
 997                                  (,op result (fast-&rest-nth i integers)))
 998                          (i 1 (1+ i)))
 999                         ((>= i (length integers))
1000                          result)
1001                       (declare (integer result)))
1002                     ,init))))
1003   (def logior 0 "Return the bit-wise or of its arguments. Args must be integers.")
1004   (def logxor 0 "Return the bit-wise exclusive or of its arguments. Args must be integers.")
1005   (def logand -1 "Return the bit-wise and of its arguments. Args must be integers.")
1006   (def logeqv -1 "Return the bit-wise equivalence of its arguments. Args must be integers."))
1007
1008 (defun lognot (number)
1009   #!+sb-doc
1010   "Return the bit-wise logical not of integer."
1011   (declare (explicit-check))
1012   (etypecase number
1013     (fixnum (lognot (truly-the fixnum number)))
1014     (bignum (bignum-logical-not number))))
1015
1016 (macrolet ((def (name explicit-check op big-op &optional doc)
1017              `(defun ,name (integer1 integer2)
1018                 ,@(when doc (list doc))
1019                 ,@(when explicit-check `((declare (explicit-check))))
1020                 (let ((x integer1)
1021                       (y integer2))
1022                   (number-dispatch ((x integer) (y integer))
1023                     (bignum-cross-fixnum ,op ,big-op))))))
1024   (def two-arg-and nil logand bignum-logical-and)
1025   (def two-arg-ior nil logior bignum-logical-ior)
1026   (def two-arg-xor nil logxor bignum-logical-xor)
1027   ;; BIGNUM-LOGICAL-{AND,IOR,XOR} need not return a bignum, so must
1028   ;; call the generic LOGNOT...
1029   (def two-arg-eqv nil logeqv (lambda (x y) (lognot (bignum-logical-xor x y))))
1030   (def lognand t lognand
1031        (lambda (x y) (lognot (bignum-logical-and x y)))
1032        #!+sb-doc "Complement the logical AND of INTEGER1 and INTEGER2.")
1033   (def lognor t lognor
1034        (lambda (x y) (lognot (bignum-logical-ior x y)))
1035        #!+sb-doc "Complement the logical OR of INTEGER1 and INTEGER2.")
1036   ;; ... but BIGNUM-LOGICAL-NOT on a bignum will always return a bignum
1037   (def logandc1 t logandc1
1038        (lambda (x y) (bignum-logical-and (bignum-logical-not x) y))
1039        #!+sb-doc "Bitwise AND (LOGNOT INTEGER1) with INTEGER2.")
1040   (def logandc2 t logandc2
1041        (lambda (x y) (bignum-logical-and x (bignum-logical-not y)))
1042        #!+sb-doc "Bitwise AND INTEGER1 with (LOGNOT INTEGER2).")
1043   (def logorc1 t logorc1
1044        (lambda (x y) (bignum-logical-ior (bignum-logical-not x) y))
1045        #!+sb-doc "Bitwise OR (LOGNOT INTEGER1) with INTEGER2.")
1046   (def logorc2 t logorc2
1047        (lambda (x y) (bignum-logical-ior x (bignum-logical-not y)))
1048        #!+sb-doc "Bitwise OR INTEGER1 with (LOGNOT INTEGER2)."))
1049
1050 (defun logcount (integer)
1051   #!+sb-doc
1052   "Count the number of 1 bits if INTEGER is non-negative,
1053 and the number of 0 bits if INTEGER is negative."
1054   (declare (explicit-check))
1055   (etypecase integer
1056     (fixnum
1057      (logcount (truly-the (integer 0
1058                                    #.(max sb!xc:most-positive-fixnum
1059                                           (lognot sb!xc:most-negative-fixnum)))
1060                           (if (minusp (truly-the fixnum integer))
1061                               (lognot (truly-the fixnum integer))
1062                               integer))))
1063     (bignum
1064      (bignum-logcount integer))))
1065
1066 (defun logtest (integer1 integer2)
1067   #!+sb-doc
1068   "Predicate which returns T if logand of integer1 and integer2 is not zero."
1069   (logtest integer1 integer2))
1070
1071 (defun logbitp (index integer)
1072   #!+sb-doc
1073   "Predicate returns T if bit index of integer is a 1."
1074   (number-dispatch ((index integer) (integer integer))
1075     ((fixnum fixnum) (if (< index sb!vm:n-positive-fixnum-bits)
1076                          (not (zerop (logand integer (ash 1 index))))
1077                          (minusp integer)))
1078     ((fixnum bignum) (bignum-logbitp index integer))
1079     ((bignum (foreach fixnum bignum)) (minusp integer))))
1080
1081 (defun ash (integer count)
1082   #!+sb-doc
1083   "Shifts integer left by count places preserving sign. - count shifts right."
1084   (declare (integer integer count) (explicit-check))
1085   (etypecase integer
1086     (fixnum
1087      (cond ((zerop integer)
1088             0)
1089            ((fixnump count)
1090             (let ((length (integer-length (truly-the fixnum integer)))
1091                   (count (truly-the fixnum count)))
1092               (declare (fixnum length count))
1093               (cond ((and (plusp count)
1094                           (>= (+ length count)
1095                               sb!vm:n-word-bits))
1096                      (bignum-ashift-left (make-small-bignum integer) count))
1097                     (t
1098                      (truly-the (signed-byte #.sb!vm:n-word-bits)
1099                                 (ash (truly-the fixnum integer) count))))))
1100            ((minusp count)
1101             (if (minusp integer) -1 0))
1102            (t
1103             (bignum-ashift-left (make-small-bignum integer) count))))
1104     (bignum
1105      (if (plusp count)
1106          (bignum-ashift-left integer count)
1107          (bignum-ashift-right integer (- count))))))
1108
1109 (defun integer-length (integer)
1110   #!+sb-doc
1111   "Return the number of non-sign bits in the twos-complement representation
1112   of INTEGER."
1113   (declare (explicit-check))
1114   (etypecase integer
1115     (fixnum
1116      (integer-length (truly-the fixnum integer)))
1117     (bignum
1118      (bignum-integer-length integer))))
1119 \f
1120 ;;;; BYTE, bytespecs, and related operations
1121
1122 (defun byte (size position)
1123   #!+sb-doc
1124   "Return a byte specifier which may be used by other byte functions
1125   (e.g. LDB)."
1126   (byte size position))
1127
1128 (defun byte-size (bytespec)
1129   #!+sb-doc
1130   "Return the size part of the byte specifier bytespec."
1131   (byte-size bytespec))
1132
1133 (defun byte-position (bytespec)
1134   #!+sb-doc
1135   "Return the position part of the byte specifier bytespec."
1136   (byte-position bytespec))
1137
1138 (defun ldb (bytespec integer)
1139   #!+sb-doc
1140   "Extract the specified byte from integer, and right justify result."
1141   (ldb bytespec integer))
1142
1143 (defun ldb-test (bytespec integer)
1144   #!+sb-doc
1145   "Return T if any of the specified bits in integer are 1's."
1146   (ldb-test bytespec integer))
1147
1148 (defun mask-field (bytespec integer)
1149   #!+sb-doc
1150   "Extract the specified byte from integer,  but do not right justify result."
1151   (mask-field bytespec integer))
1152
1153 (defun dpb (newbyte bytespec integer)
1154   #!+sb-doc
1155   "Return new integer with newbyte in specified position, newbyte is right justified."
1156   (dpb newbyte bytespec integer))
1157
1158 (defun deposit-field (newbyte bytespec integer)
1159   #!+sb-doc
1160   "Return new integer with newbyte in specified position, newbyte is not right justified."
1161   (deposit-field newbyte bytespec integer))
1162
1163 (defun %ldb (size posn integer)
1164   (declare (type bit-index size posn) (explicit-check))
1165   ;; The naive algorithm is horrible in the general case.
1166   ;; Consider (LDB (BYTE 1 2) (SOME-GIANT-BIGNUM)) which has to shift the
1167   ;; input rightward 2 bits, consing a new bignum just to read 1 bit.
1168   (if (and (<= 0 size sb!vm:n-positive-fixnum-bits)
1169            (typep integer 'bignum))
1170       (sb!bignum::ldb-bignum=>fixnum size posn integer)
1171       (logand (ash integer (- posn))
1172               (1- (ash 1 size)))))
1173
1174 (defun %mask-field (size posn integer)
1175   (declare (type bit-index size posn) (explicit-check))
1176   (logand integer (ash (1- (ash 1 size)) posn)))
1177
1178 (defun %dpb (newbyte size posn integer)
1179   (declare (type bit-index size posn) (explicit-check))
1180   (let ((mask (1- (ash 1 size))))
1181     (logior (logand integer (lognot (ash mask posn)))
1182             (ash (logand newbyte mask) posn))))
1183
1184 (defun %deposit-field (newbyte size posn integer)
1185   (declare (type bit-index size posn) (explicit-check))
1186   (let ((mask (ash (ldb (byte size 0) -1) posn)))
1187     (logior (logand newbyte mask)
1188             (logand integer (lognot mask)))))
1189
1190 (defun sb!c::mask-signed-field (size integer)
1191   #!+sb-doc
1192   "Extract SIZE lower bits from INTEGER, considering them as a
1193 2-complement SIZE-bits representation of a signed integer."
1194   (macrolet ((msf (size integer)
1195                `(if (logbitp (1- ,size) ,integer)
1196                     (dpb ,integer (byte (1- ,size) 0) -1)
1197                     (ldb (byte (1- ,size) 0) ,integer))))
1198     (typecase size
1199       ((eql 0) 0)
1200       ((integer 1 #.sb!vm:n-fixnum-bits)
1201        (number-dispatch ((integer integer))
1202          ((fixnum) (msf size integer))
1203          ((bignum) (let ((fix (sb!c::mask-signed-field #.sb!vm:n-fixnum-bits (%bignum-ref integer 0))))
1204                      (if (= size #.sb!vm:n-fixnum-bits)
1205                          fix
1206                          (msf size fix))))))
1207       ((integer (#.sb!vm:n-fixnum-bits) #.sb!vm:n-word-bits)
1208        (number-dispatch ((integer integer))
1209          ((fixnum) integer)
1210          ((bignum) (let ((word (sb!c::mask-signed-field #.sb!vm:n-word-bits (%bignum-ref integer 0))))
1211                      (if (= size #.sb!vm:n-word-bits)
1212                          word
1213                          (msf size word))))))
1214       ((unsigned-byte) (msf size integer)))))
1215 \f
1216 ;;;; BOOLE
1217
1218 ;;; The boole function dispaches to any logic operation depending on
1219 ;;;     the value of a variable. Presently, legal selector values are [0..15].
1220 ;;;     boole is open coded for calls with a constant selector. or with calls
1221 ;;;     using any of the constants declared below.
1222
1223 (defconstant boole-clr 0
1224   #!+sb-doc
1225   "Boole function op, makes BOOLE return 0.")
1226
1227 (defconstant boole-set 1
1228   #!+sb-doc
1229   "Boole function op, makes BOOLE return -1.")
1230
1231 (defconstant boole-1   2
1232   #!+sb-doc
1233   "Boole function op, makes BOOLE return integer1.")
1234
1235 (defconstant boole-2   3
1236   #!+sb-doc
1237   "Boole function op, makes BOOLE return integer2.")
1238
1239 (defconstant boole-c1  4
1240   #!+sb-doc
1241   "Boole function op, makes BOOLE return complement of integer1.")
1242
1243 (defconstant boole-c2  5
1244   #!+sb-doc
1245   "Boole function op, makes BOOLE return complement of integer2.")
1246
1247 (defconstant boole-and 6
1248   #!+sb-doc
1249   "Boole function op, makes BOOLE return logand of integer1 and integer2.")
1250
1251 (defconstant boole-ior 7
1252   #!+sb-doc
1253   "Boole function op, makes BOOLE return logior of integer1 and integer2.")
1254
1255 (defconstant boole-xor 8
1256   #!+sb-doc
1257   "Boole function op, makes BOOLE return logxor of integer1 and integer2.")
1258
1259 (defconstant boole-eqv 9
1260   #!+sb-doc
1261   "Boole function op, makes BOOLE return logeqv of integer1 and integer2.")
1262
1263 (defconstant boole-nand  10
1264   #!+sb-doc
1265   "Boole function op, makes BOOLE return log nand of integer1 and integer2.")
1266
1267 (defconstant boole-nor   11
1268   #!+sb-doc
1269   "Boole function op, makes BOOLE return lognor of integer1 and integer2.")
1270
1271 (defconstant boole-andc1 12
1272   #!+sb-doc
1273   "Boole function op, makes BOOLE return logandc1 of integer1 and integer2.")
1274
1275 (defconstant boole-andc2 13
1276   #!+sb-doc
1277   "Boole function op, makes BOOLE return logandc2 of integer1 and integer2.")
1278
1279 (defconstant boole-orc1  14
1280   #!+sb-doc
1281   "Boole function op, makes BOOLE return logorc1 of integer1 and integer2.")
1282
1283 (defconstant boole-orc2  15
1284   #!+sb-doc
1285   "Boole function op, makes BOOLE return logorc2 of integer1 and integer2.")
1286
1287 (defun boole (op integer1 integer2)
1288   #!+sb-doc
1289   "Bit-wise boolean function on two integers. Function chosen by OP:
1290         0       BOOLE-CLR
1291         1       BOOLE-SET
1292         2       BOOLE-1
1293         3       BOOLE-2
1294         4       BOOLE-C1
1295         5       BOOLE-C2
1296         6       BOOLE-AND
1297         7       BOOLE-IOR
1298         8       BOOLE-XOR
1299         9       BOOLE-EQV
1300         10      BOOLE-NAND
1301         11      BOOLE-NOR
1302         12      BOOLE-ANDC1
1303         13      BOOLE-ANDC2
1304         14      BOOLE-ORC1
1305         15      BOOLE-ORC2"
1306   (case op
1307     (0 (boole 0 integer1 integer2))
1308     (1 (boole 1 integer1 integer2))
1309     (2 (boole 2 integer1 integer2))
1310     (3 (boole 3 integer1 integer2))
1311     (4 (boole 4 integer1 integer2))
1312     (5 (boole 5 integer1 integer2))
1313     (6 (boole 6 integer1 integer2))
1314     (7 (boole 7 integer1 integer2))
1315     (8 (boole 8 integer1 integer2))
1316     (9 (boole 9 integer1 integer2))
1317     (10 (boole 10 integer1 integer2))
1318     (11 (boole 11 integer1 integer2))
1319     (12 (boole 12 integer1 integer2))
1320     (13 (boole 13 integer1 integer2))
1321     (14 (boole 14 integer1 integer2))
1322     (15 (boole 15 integer1 integer2))
1323     (t (error 'type-error :datum op :expected-type '(mod 16)))))
1324 \f
1325 ;;;; GCD and LCM
1326
1327 (defun gcd (&rest integers)
1328   #!+sb-doc
1329   "Return the greatest common divisor of the arguments, which must be
1330   integers. GCD with no arguments is defined to be 0."
1331   (declare (explicit-check))
1332   (case (length integers)
1333     (0 0)
1334     (1 (abs (the integer (fast-&rest-nth 0 integers))))
1335     (otherwise
1336      (do ((result (fast-&rest-nth 0 integers)
1337                   (gcd result (the integer (fast-&rest-nth i integers))))
1338           (i 1 (1+ i)))
1339          ((>= i (length integers))
1340           result)
1341        (declare (integer result))))))
1342
1343 (defun lcm (&rest integers)
1344   #!+sb-doc
1345   "Return the least common multiple of one or more integers. LCM of no
1346   arguments is defined to be 1."
1347   (declare (explicit-check))
1348   (case (length integers)
1349     (0 1)
1350     (1 (abs (the integer (fast-&rest-nth 0 integers))))
1351     (otherwise
1352      (do ((result (fast-&rest-nth 0 integers)
1353                   (lcm result (the integer (fast-&rest-nth i integers))))
1354           (i 1 (1+ i)))
1355          ((>= i (length integers))
1356           result)
1357        (declare (integer result))))))
1358
1359 (defun two-arg-lcm (n m)
1360   (declare (integer n m))
1361   (if (or (zerop n) (zerop m))
1362       0
1363       ;; KLUDGE: I'm going to assume that it was written this way
1364       ;; originally for a reason.  However, this is a somewhat
1365       ;; complicated way of writing the algorithm in the CLHS page for
1366       ;; LCM, and I don't know why.  To be investigated.  -- CSR,
1367       ;; 2003-09-11
1368       ;;
1369       ;;    It seems to me that this is written this way to avoid
1370       ;;    unnecessary bignumification of intermediate results.
1371       ;;        -- TCR, 2008-03-05
1372       (let ((m (abs m))
1373             (n (abs n)))
1374         (multiple-value-bind (max min)
1375             (if (> m n)
1376                 (values m n)
1377                 (values n m))
1378           (* (truncate max (gcd n m)) min)))))
1379
1380 ;;; Do the GCD of two integer arguments. With fixnum arguments, we use the
1381 ;;; binary GCD algorithm from Knuth's seminumerical algorithms (slightly
1382 ;;; structurified), otherwise we call BIGNUM-GCD. We pick off the special case
1383 ;;; of 0 before the dispatch so that the bignum code doesn't have to worry
1384 ;;; about "small bignum" zeros.
1385 (defun two-arg-gcd (u v)
1386   (cond ((eql u 0) (abs v))
1387         ((eql v 0) (abs u))
1388         (t
1389          (number-dispatch ((u integer) (v integer))
1390            ((fixnum fixnum)
1391             (locally
1392                 (declare (optimize (speed 3) (safety 0)))
1393               (do ((k 0 (1+ k))
1394                    (u (abs u) (ash u -1))
1395                    (v (abs v) (ash v -1)))
1396                   ((oddp (logior u v))
1397                    (do ((temp (if (oddp u) (- v) (ash u -1))
1398                               (ash temp -1)))
1399                        (nil)
1400                      (declare (fixnum temp))
1401                      (when (oddp temp)
1402                        (if (plusp temp)
1403                            (setq u temp)
1404                            (setq v (- temp)))
1405                        (setq temp (- u v))
1406                        (when (zerop temp)
1407                          (let ((res (ash u k)))
1408                            (declare (type sb!vm:signed-word res)
1409                                     (optimize (inhibit-warnings 3)))
1410                            (return res))))))
1411                 (declare (type (mod #.sb!vm:n-word-bits) k)
1412                          (type sb!vm:signed-word u v)))))
1413            ((bignum bignum)
1414             (bignum-gcd u v))
1415            ((bignum fixnum)
1416             (bignum-gcd u (make-small-bignum v)))
1417            ((fixnum bignum)
1418             (bignum-gcd (make-small-bignum u) v))))))
1419 \f
1420 ;;; from Robert Smith; changed not to cons unnecessarily, and tuned for
1421 ;;; faster operation on fixnum inputs by compiling the central recursive
1422 ;;; algorithm twice, once using generic and once fixnum arithmetic, and
1423 ;;; dispatching on function entry into the applicable part. For maximum
1424 ;;; speed, the fixnum part recurs into itself, thereby avoiding further
1425 ;;; type dispatching. This pattern is not supported by NUMBER-DISPATCH
1426 ;;; thus some special-purpose macrology is needed.
1427 (defun isqrt (n)
1428   #!+sb-doc
1429   "Return the greatest integer less than or equal to the square root of N."
1430   (declare (type unsigned-byte n) (explicit-check))
1431   (macrolet
1432       ((isqrt-recursion (arg recurse fixnum-p)
1433          ;; Expands into code for the recursive step of the ISQRT
1434          ;; calculation. ARG is the input variable and RECURSE the name
1435          ;; of the function to recur into. If FIXNUM-P is true, some
1436          ;; type declarations are added that, together with ARG being
1437          ;; declared as a fixnum outside of here, make the resulting code
1438          ;; compile into fixnum-specialized code without any calls to
1439          ;; generic arithmetic. Else, the code works for bignums, too.
1440          ;; The input must be at least 16 to ensure that RECURSE is called
1441          ;; with a strictly smaller number and that the result is correct
1442          ;; (provided that RECURSE correctly implements ISQRT, itself).
1443          `(macrolet ((if-fixnum-p-truly-the (type expr)
1444                        ,@(if fixnum-p
1445                              '(`(truly-the ,type ,expr))
1446                              '((declare (ignore type))
1447                                expr))))
1448             (let* ((fourth-size (ash (1- (integer-length ,arg)) -2))
1449                    (significant-half (ash ,arg (- (ash fourth-size 1))))
1450                    (significant-half-isqrt
1451                     (if-fixnum-p-truly-the
1452                      (integer 1 #.(isqrt sb!xc:most-positive-fixnum))
1453                      (,recurse significant-half)))
1454                    (zeroth-iteration (ash significant-half-isqrt
1455                                           fourth-size)))
1456               (multiple-value-bind (quot rem)
1457                   (floor ,arg zeroth-iteration)
1458                 (let ((first-iteration (ash (+ zeroth-iteration quot) -1)))
1459                   (cond ((oddp quot)
1460                          first-iteration)
1461                         ((> (if-fixnum-p-truly-the
1462                              fixnum
1463                              (expt (- first-iteration zeroth-iteration) 2))
1464                             rem)
1465                          (1- first-iteration))
1466                         (t
1467                          first-iteration))))))))
1468     (typecase n
1469       (fixnum (labels ((fixnum-isqrt (n)
1470                          (declare (type fixnum n))
1471                          (cond ((> n 24)
1472                                 (isqrt-recursion n fixnum-isqrt t))
1473                                ((> n 15) 4)
1474                                ((> n  8) 3)
1475                                ((> n  3) 2)
1476                                ((> n  0) 1)
1477                                ((= n  0) 0))))
1478                 (fixnum-isqrt n)))
1479       (bignum (isqrt-recursion n isqrt nil)))))
1480 \f
1481 ;;;; miscellaneous number predicates
1482
1483 (macrolet ((def (name doc)
1484              (declare (ignorable doc))
1485              `(defun ,name (number) #!+sb-doc ,doc
1486                 (declare (explicit-check))
1487                 (,name number))))
1488   (def zerop "Is this number zero?")
1489   (def plusp "Is this real number strictly positive?")
1490   (def minusp "Is this real number strictly negative?")
1491   (def oddp "Is this integer odd?")
1492   (def evenp "Is this integer even?"))
1493 \f
1494 ;;;; modular functions
1495 #.
1496 (collect ((forms))
1497   (flet ((unsigned-definition (name lambda-list width)
1498            (let ((pattern (1- (ash 1 width))))
1499              `(defun ,name ,(copy-list lambda-list)
1500                (flet ((prepare-argument (x)
1501                         (declare (integer x))
1502                         (etypecase x
1503                           ((unsigned-byte ,width) x)
1504                           (fixnum (logand x ,pattern))
1505                           (bignum (logand x ,pattern)))))
1506                  (,name ,@(loop for arg in lambda-list
1507                                 collect `(prepare-argument ,arg)))))))
1508          (signed-definition (name lambda-list width)
1509            `(defun ,name ,(copy-list lambda-list)
1510               (flet ((prepare-argument (x)
1511                        (declare (integer x))
1512                        (etypecase x
1513                          ((signed-byte ,width) x)
1514                          (fixnum (sb!c::mask-signed-field ,width x))
1515                          (bignum (sb!c::mask-signed-field ,width x)))))
1516                 (,name ,@(loop for arg in lambda-list
1517                                collect `(prepare-argument ,arg)))))))
1518     (flet ((do-mfuns (class)
1519              (loop for infos being each hash-value of (sb!c::modular-class-funs class)
1520                    ;; FIXME: We need to process only "toplevel" functions
1521                    when (listp infos)
1522                    do (loop for info in infos
1523                             for name = (sb!c::modular-fun-info-name info)
1524                             and width = (sb!c::modular-fun-info-width info)
1525                             and signedp = (sb!c::modular-fun-info-signedp info)
1526                             and lambda-list = (sb!c::modular-fun-info-lambda-list info)
1527                             if signedp
1528                             do (forms (signed-definition name lambda-list width))
1529                             else
1530                             do (forms (unsigned-definition name lambda-list width))))))
1531       (do-mfuns sb!c::*untagged-unsigned-modular-class*)
1532       (do-mfuns sb!c::*untagged-signed-modular-class*)
1533       (do-mfuns sb!c::*tagged-modular-class*)))
1534   `(progn ,@(sort (forms) #'string< :key #'cadr)))
1535
1536 ;;; KLUDGE: these out-of-line definitions can't use the modular
1537 ;;; arithmetic, as that is only (currently) defined for constant
1538 ;;; shifts.  See also the comment in (LOGAND OPTIMIZER) for more
1539 ;;; discussion of this hack.  -- CSR, 2003-10-09
1540 #!-64-bit-registers
1541 (defun sb!vm::ash-left-mod32 (integer amount)
1542   (etypecase integer
1543     ((unsigned-byte 32) (ldb (byte 32 0) (ash integer amount)))
1544     (fixnum (ldb (byte 32 0) (ash (logand integer #xffffffff) amount)))
1545     (bignum (ldb (byte 32 0) (ash (logand integer #xffffffff) amount)))))
1546 #!+64-bit-registers
1547 (defun sb!vm::ash-left-mod64 (integer amount)
1548   (etypecase integer
1549     ((unsigned-byte 64) (ldb (byte 64 0) (ash integer amount)))
1550     (fixnum (ldb (byte 64 0) (ash (logand integer #xffffffffffffffff) amount)))
1551     (bignum (ldb (byte 64 0)
1552                  (ash (logand integer #xffffffffffffffff) amount)))))
1553
1554 #!+(or x86 x86-64 arm arm64)
1555 (defun sb!vm::ash-left-modfx (integer amount)
1556   (let ((fixnum-width (- sb!vm:n-word-bits sb!vm:n-fixnum-tag-bits)))
1557     (etypecase integer
1558       (fixnum (sb!c::mask-signed-field fixnum-width (ash integer amount)))
1559       (integer (sb!c::mask-signed-field fixnum-width (ash (sb!c::mask-signed-field fixnum-width integer) amount))))))