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   (flet ((%%make-complex (realpart imagpart)
 266            (cond #!+long-float
 267                  ((and (typep realpart 'long-float)
 268                        (typep imagpart 'long-float))
 269                   (truly-the (complex long-float)
 270                              (complex realpart imagpart)))
 271                  ((and (typep realpart 'double-float)
 272                        (typep imagpart 'double-float))
 273                   (truly-the (complex double-float)
 274                              (complex realpart imagpart)))
 275                  ((and (typep realpart 'single-float)
 276                        (typep imagpart 'single-float))
 277                   (truly-the (complex single-float)
 278                              (complex realpart imagpart)))
 279                  (t
 280                   (%make-complex realpart imagpart)))))
 281   (number-dispatch ((realpart real) (imagpart real))
 282     ((rational rational)
 283      (canonical-complex realpart imagpart))
 284     (float-contagion %%make-complex realpart imagpart (rational)))))
 285
 286 (defun realpart (number)
 287   #!+sb-doc
 288   "Extract the real part of a number."
 289   (etypecase number
 290     #!+long-float
 291     ((complex long-float)
 292      (truly-the long-float (realpart number)))
 293     ((complex double-float)
 294      (truly-the double-float (realpart number)))
 295     ((complex single-float)
 296      (truly-the single-float (realpart number)))
 297     ((complex rational)
 298      (sb!kernel:%realpart number))
 299     (number
 300      number)))
 301
 302 (defun imagpart (number)
 303   #!+sb-doc
 304   "Extract the imaginary part of a number."
 305   (etypecase number
 306     #!+long-float
 307     ((complex long-float)
 308      (truly-the long-float (imagpart number)))
 309     ((complex double-float)
 310      (truly-the double-float (imagpart number)))
 311     ((complex single-float)
 312      (truly-the single-float (imagpart number)))
 313     ((complex rational)
 314      (sb!kernel:%imagpart number))
 315     (float
 316      (* 0 number))
 317     (number
 318      0)))
 319
 320 (defun conjugate (number)
 321   #!+sb-doc
 322   "Return the complex conjugate of NUMBER. For non-complex numbers, this is
 323   an identity."
 324   (declare (type number number))
 325   (if (complexp number)
 326       (complex (realpart number) (- (imagpart number)))
 327       number))
 328
 329 (defun signum (number)
 330   #!+sb-doc
 331   "If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))."
 332   (if (zerop number)
 333       number
 334       (if (rationalp number)
 335           (if (plusp number) 1 -1)
 336           (/ number (abs number)))))
 337 \f
 338 ;;;; ratios
 339
 340 (defun numerator (number)
 341   #!+sb-doc
 342   "Return the numerator of NUMBER, which must be rational."
 343   (numerator number))
 344
 345 (defun denominator (number)
 346   #!+sb-doc
 347   "Return the denominator of NUMBER, which must be rational."
 348   (denominator number))
 349 \f
 350 ;;;; arithmetic operations
 351
 352 (macrolet ((define-arith (op init doc)
 353              #!-sb-doc (declare (ignore doc))
 354              `(defun ,op (&rest args)
 355                 #!+sb-doc ,doc
 356                 (declare (truly-dynamic-extent args))
 357                 (if (null args) ,init
 358                     (do ((args (cdr args) (cdr args))
 359                          (result (car args) (,op result (car args))))
 360                         ((null args) result)
 361                       ;; to signal TYPE-ERROR when exactly 1 arg of wrong type:
 362                       (declare (type number result)))))))
 363   (define-arith + 0
 364     "Return the sum of its arguments. With no args, returns 0.")
 365   (define-arith * 1
 366     "Return the product of its arguments. With no args, returns 1."))
 367
 368 (defun - (number &rest more-numbers)
 369   #!+sb-doc
 370   "Subtract the second and all subsequent arguments from the first;
 371   or with one argument, negate the first argument."
 372   (declare (truly-dynamic-extent more-numbers))
 373   (if more-numbers
 374       (do ((nlist more-numbers (cdr nlist))
 375            (result number))
 376           ((atom nlist) result)
 377          (declare (list nlist))
 378          (setq result (- result (car nlist))))
 379       (- number)))
 380
 381 (defun / (number &rest more-numbers)
 382   #!+sb-doc
 383   "Divide the first argument by each of the following arguments, in turn.
 384   With one argument, return reciprocal."
 385   (declare (truly-dynamic-extent more-numbers))
 386   (if more-numbers
 387       (do ((nlist more-numbers (cdr nlist))
 388            (result number))
 389           ((atom nlist) result)
 390          (declare (list nlist))
 391          (setq result (/ result (car nlist))))
 392       (/ number)))
 393
 394 (defun 1+ (number)
 395   #!+sb-doc
 396   "Return NUMBER + 1."
 397   (1+ number))
 398
 399 (defun 1- (number)
 400   #!+sb-doc
 401   "Return NUMBER - 1."
 402   (1- number))
 403
 404 (eval-when (:compile-toplevel)
 405
 406 (sb!xc:defmacro two-arg-+/- (name op big-op)
 407   `(defun ,name (x y)
 408      (number-dispatch ((x number) (y number))
 409        (bignum-cross-fixnum ,op ,big-op)
 410        (float-contagion ,op x y)
 411
 412        ((complex complex)
 413         (canonical-complex (,op (realpart x) (realpart y))
 414                            (,op (imagpart x) (imagpart y))))
 415        (((foreach bignum fixnum ratio single-float double-float
 416                   #!+long-float long-float) complex)
 417         (complex (,op x (realpart y)) (,op 0 (imagpart y))))
 418        ((complex (or rational float))
 419         (complex (,op (realpart x) y) (,op (imagpart x) 0)))
 420
 421        (((foreach fixnum bignum) ratio)
 422         (let* ((dy (denominator y))
 423                (n (,op (* x dy) (numerator y))))
 424           (%make-ratio n dy)))
 425        ((ratio integer)
 426         (let* ((dx (denominator x))
 427                (n (,op (numerator x) (* y dx))))
 428           (%make-ratio n dx)))
 429        ((ratio ratio)
 430         (let* ((nx (numerator x))
 431                (dx (denominator x))
 432                (ny (numerator y))
 433                (dy (denominator y))
 434                (g1 (gcd dx dy)))
 435           (if (eql g1 1)
 436               (%make-ratio (,op (* nx dy) (* dx ny)) (* dx dy))
 437               (let* ((t1 (,op (* nx (truncate dy g1)) (* (truncate dx g1) ny)))
 438                      (g2 (gcd t1 g1))
 439                      (t2 (truncate dx g1)))
 440                 (cond ((eql t1 0) 0)
 441                       ((eql g2 1)
 442                        (%make-ratio t1 (* t2 dy)))
 443                       (t (let* ((nn (truncate t1 g2))
 444                                 (t3 (truncate dy g2))
 445                                 (nd (if (eql t2 1) t3 (* t2 t3))))
 446                            (if (eql nd 1) nn (%make-ratio nn nd))))))))))))
 447
 448 ) ; EVAL-WHEN
 449
 450 (two-arg-+/- two-arg-+ + add-bignums)
 451 (two-arg-+/- two-arg-- - subtract-bignum)
 452
 453 (defun two-arg-* (x y)
 454   (flet ((integer*ratio (x y)
 455            (if (eql x 0) 0
 456                (let* ((ny (numerator y))
 457                       (dy (denominator y))
 458                       (gcd (gcd x dy)))
 459                  (if (eql gcd 1)
 460                      (%make-ratio (* x ny) dy)
 461                      (let ((nn (* (truncate x gcd) ny))
 462                            (nd (truncate dy gcd)))
 463                        (if (eql nd 1)
 464                            nn
 465                            (%make-ratio nn nd)))))))
 466          (complex*real (x y)
 467            (canonical-complex (* (realpart x) y) (* (imagpart x) y))))
 468     (number-dispatch ((x number) (y number))
 469       (float-contagion * x y)
 470
 471       ((fixnum fixnum) (multiply-fixnums x y))
 472       ((bignum fixnum) (multiply-bignum-and-fixnum x y))
 473       ((fixnum bignum) (multiply-bignum-and-fixnum y x))
 474       ((bignum bignum) (multiply-bignums x y))
 475
 476       ((complex complex)
 477        (let* ((rx (realpart x))
 478               (ix (imagpart x))
 479               (ry (realpart y))
 480               (iy (imagpart y)))
 481          (canonical-complex (- (* rx ry) (* ix iy)) (+ (* rx iy) (* ix ry)))))
 482       (((foreach bignum fixnum ratio single-float double-float
 483                  #!+long-float long-float)
 484         complex)
 485        (complex*real y x))
 486       ((complex (or rational float))
 487        (complex*real x y))
 488
 489       (((foreach bignum fixnum) ratio) (integer*ratio x y))
 490       ((ratio integer) (integer*ratio y x))
 491       ((ratio ratio)
 492        (let* ((nx (numerator x))
 493               (dx (denominator x))
 494               (ny (numerator y))
 495               (dy (denominator y))
 496               (g1 (gcd nx dy))
 497               (g2 (gcd dx ny)))
 498          (build-ratio (* (maybe-truncate nx g1)
 499                          (maybe-truncate ny g2))
 500                       (* (maybe-truncate dx g2)
 501                          (maybe-truncate dy g1))))))))
 502
 503 ;;; Divide two integers, producing a canonical rational. If a fixnum,
 504 ;;; we see whether they divide evenly before trying the GCD. In the
 505 ;;; bignum case, we don't bother, since bignum division is expensive,
 506 ;;; and the test is not very likely to succeed.
 507 (defun integer-/-integer (x y)
 508   (if (and (typep x 'fixnum) (typep y 'fixnum))
 509       (multiple-value-bind (quo rem) (truncate x y)
 510         (if (zerop rem)
 511             quo
 512             (let ((gcd (gcd x y)))
 513               (declare (fixnum gcd))
 514               (if (eql gcd 1)
 515                   (build-ratio x y)
 516                   (build-ratio (truncate x gcd) (truncate y gcd))))))
 517       (let ((gcd (gcd x y)))
 518         (if (eql gcd 1)
 519             (build-ratio x y)
 520             (build-ratio (truncate x gcd) (truncate y gcd))))))
 521
 522 (defun two-arg-/ (x y)
 523   (number-dispatch ((x number) (y number))
 524     (float-contagion / x y (ratio integer))
 525
 526     ((complex complex)
 527      (let* ((rx (realpart x))
 528             (ix (imagpart x))
 529             (ry (realpart y))
 530             (iy (imagpart y)))
 531        (if (> (abs ry) (abs iy))
 532            (let* ((r (/ iy ry))
 533                   (dn (* ry (+ 1 (* r r)))))
 534              (canonical-complex (/ (+ rx (* ix r)) dn)
 535                                 (/ (- ix (* rx r)) dn)))
 536            (let* ((r (/ ry iy))
 537                   (dn (* iy (+ 1 (* r r)))))
 538              (canonical-complex (/ (+ (* rx r) ix) dn)
 539                                 (/ (- (* ix r) rx) dn))))))
 540     (((foreach integer ratio single-float double-float) complex)
 541      (let* ((ry (realpart y))
 542             (iy (imagpart y)))
 543        (if (> (abs ry) (abs iy))
 544            (let* ((r (/ iy ry))
 545                   (dn (* ry (+ 1 (* r r)))))
 546              (canonical-complex (/ x dn)
 547                                 (/ (- (* x r)) dn)))
 548            (let* ((r (/ ry iy))
 549                   (dn (* iy (+ 1 (* r r)))))
 550              (canonical-complex (/ (* x r) dn)
 551                                 (/ (- x) dn))))))
 552     ((complex (or rational float))
 553      (canonical-complex (/ (realpart x) y)
 554                         (/ (imagpart x) y)))
 555
 556     ((ratio ratio)
 557      (let* ((nx (numerator x))
 558             (dx (denominator x))
 559             (ny (numerator y))
 560             (dy (denominator y))
 561             (g1 (gcd nx ny))
 562             (g2 (gcd dx dy)))
 563        (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2))
 564                     (* (maybe-truncate dx g2) (maybe-truncate ny g1)))))
 565
 566     ((integer integer)
 567      (integer-/-integer x y))
 568
 569     ((integer ratio)
 570      (if (zerop x)
 571          0
 572          (let* ((ny (numerator y))
 573                 (dy (denominator y))
 574                 (gcd (gcd x ny)))
 575            (build-ratio (* (maybe-truncate x gcd) dy)
 576                         (maybe-truncate ny gcd)))))
 577
 578     ((ratio integer)
 579      (let* ((nx (numerator x))
 580             (gcd (gcd nx y)))
 581        (build-ratio (maybe-truncate nx gcd)
 582                     (* (maybe-truncate y gcd) (denominator x)))))))
 583
 584 (defun %negate (n)
 585   (number-dispatch ((n number))
 586     (((foreach fixnum single-float double-float #!+long-float long-float))
 587      (%negate n))
 588     ((bignum)
 589      (negate-bignum n))
 590     ((ratio)
 591      (%make-ratio (- (numerator n)) (denominator n)))
 592     ((complex)
 593      (complex (- (realpart n)) (- (imagpart n))))))
 594 \f
 595 ;;;; TRUNCATE and friends
 596
 597 (defun truncate (number &optional (divisor 1))
 598   #!+sb-doc
 599   "Return number (or number/divisor) as an integer, rounded toward 0.
 600   The second returned value is the remainder."
 601   (macrolet ((truncate-float (rtype)
 602                `(let* ((float-div (coerce divisor ',rtype))
 603                        (res (%unary-truncate (/ number float-div))))
 604                   (values res
 605                           (- number
 606                              (* (coerce res ',rtype) float-div))))))
 607     (number-dispatch ((number real) (divisor real))
 608       ((fixnum fixnum) (truncate number divisor))
 609       (((foreach fixnum bignum) ratio)
 610        (let ((q (truncate (* number (denominator divisor))
 611                           (numerator divisor))))
 612          (values q (- number (* q divisor)))))
 613       ((fixnum bignum)
 614        (bignum-truncate (make-small-bignum number) divisor))
 615       ((ratio (or float rational))
 616        (let ((q (truncate (numerator number)
 617                           (* (denominator number) divisor))))
 618          (values q (- number (* q divisor)))))
 619       ((bignum fixnum)
 620        (bignum-truncate number (make-small-bignum divisor)))
 621       ((bignum bignum)
 622        (bignum-truncate number divisor))
 623
 624       (((foreach single-float double-float #!+long-float long-float)
 625         (or rational single-float))
 626        (if (eql divisor 1)
 627            (let ((res (%unary-truncate number)))
 628              (values res (- number (coerce res '(dispatch-type number)))))
 629            (truncate-float (dispatch-type number))))
 630       #!+long-float
 631       ((long-float (or single-float double-float long-float))
 632        (truncate-float long-float))
 633       #!+long-float
 634       (((foreach double-float single-float) long-float)
 635        (truncate-float long-float))
 636       ((double-float (or single-float double-float))
 637        (truncate-float double-float))
 638       ((single-float double-float)
 639        (truncate-float double-float))
 640       (((foreach fixnum bignum ratio)
 641         (foreach single-float double-float #!+long-float long-float))
 642        (truncate-float (dispatch-type divisor))))))
 643
 644 ;; Only inline when no VOP exists
 645 #!-multiply-high-vops (declaim (inline %multiply-high))
 646 (defun %multiply-high (x y)
 647   (declare (type word x y))
 648   #!-multiply-high-vops
 649   (values (sb!bignum:%multiply x y))
 650   #!+multiply-high-vops
 651   (%multiply-high x y))
 652
 653 ;;; Declare these guys inline to let them get optimized a little.
 654 ;;; ROUND and FROUND are not declared inline since they seem too
 655 ;;; obscure and too big to inline-expand by default. Also, this gives
 656 ;;; the compiler a chance to pick off the unary float case.
 657 ;;;
 658 ;;; CEILING and FLOOR are implemented in terms of %CEILING and %FLOOR
 659 ;;; if no better transform can be found: they aren't inline directly,
 660 ;;; since we want to try a transform specific to them before letting
 661 ;;; the transform for TRUNCATE pick up the slack.
 662 #!-sb-fluid (declaim (inline rem mod fceiling ffloor ftruncate %floor %ceiling))
 663 (defun %floor (number divisor)
 664   ;; If the numbers do not divide exactly and the result of
 665   ;; (/ NUMBER DIVISOR) would be negative then decrement the quotient
 666   ;; and augment the remainder by the divisor.
 667   (multiple-value-bind (tru rem) (truncate number divisor)
 668     (if (and (not (zerop rem))
 669              (if (minusp divisor)
 670                  (plusp number)
 671                  (minusp number)))
 672         (values (1- tru) (+ rem divisor))
 673         (values tru rem))))
 674
 675 (defun floor (number &optional (divisor 1))
 676   #!+sb-doc
 677   "Return the greatest integer not greater than number, or number/divisor.
 678   The second returned value is (mod number divisor)."
 679   (%floor number divisor))
 680
 681 (defun %ceiling (number divisor)
 682   ;; If the numbers do not divide exactly and the result of
 683   ;; (/ NUMBER DIVISOR) would be positive then increment the quotient
 684   ;; and decrement the remainder by the divisor.
 685   (multiple-value-bind (tru rem) (truncate number divisor)
 686     (if (and (not (zerop rem))
 687              (if (minusp divisor)
 688                  (minusp number)
 689                  (plusp number)))
 690         (values (+ tru 1) (- rem divisor))
 691         (values tru rem))))
 692
 693 (defun ceiling (number &optional (divisor 1))
 694   #!+sb-doc
 695   "Return the smallest integer not less than number, or number/divisor.
 696   The second returned value is the remainder."
 697   (%ceiling number divisor))
 698
 699 (defun round (number &optional (divisor 1))
 700   #!+sb-doc
 701   "Rounds number (or number/divisor) to nearest integer.
 702   The second returned value is the remainder."
 703   (if (eql divisor 1)
 704       (round number)
 705       (multiple-value-bind (tru rem) (truncate number divisor)
 706         (if (zerop rem)
 707             (values tru rem)
 708             (let ((thresh (/ (abs divisor) 2)))
 709               (cond ((or (> rem thresh)
 710                          (and (= rem thresh) (oddp tru)))
 711                      (if (minusp divisor)
 712                          (values (- tru 1) (+ rem divisor))
 713                          (values (+ tru 1) (- rem divisor))))
 714                     ((let ((-thresh (- thresh)))
 715                        (or (< rem -thresh)
 716                            (and (= rem -thresh) (oddp tru))))
 717                      (if (minusp divisor)
 718                          (values (+ tru 1) (- rem divisor))
 719                          (values (- tru 1) (+ rem divisor))))
 720                     (t (values tru rem))))))))
 721
 722 (defun rem (number divisor)
 723   #!+sb-doc
 724   "Return second result of TRUNCATE."
 725   (multiple-value-bind (tru rem) (truncate number divisor)
 726     (declare (ignore tru))
 727     rem))
 728
 729 (defun mod (number divisor)
 730   #!+sb-doc
 731   "Return second result of FLOOR."
 732   (let ((rem (rem number divisor)))
 733     (if (and (not (zerop rem))
 734              (if (minusp divisor)
 735                  (plusp number)
 736                  (minusp number)))
 737         (+ rem divisor)
 738         rem)))
 739
 740 (defmacro !define-float-rounding-function (name op doc)
 741   `(defun ,name (number &optional (divisor 1))
 742     ,doc
 743     (multiple-value-bind (res rem) (,op number divisor)
 744       (values (float res (if (floatp rem) rem 1.0)) rem))))
 745
 746 (defun ftruncate (number &optional (divisor 1))
 747   #!+sb-doc
 748   "Same as TRUNCATE, but returns first value as a float."
 749   (macrolet ((ftruncate-float (rtype)
 750                `(let* ((float-div (coerce divisor ',rtype))
 751                        (res (%unary-ftruncate (/ number float-div))))
 752                   (values res
 753                           (- number
 754                              (* (coerce res ',rtype) float-div))))))
 755     (number-dispatch ((number real) (divisor real))
 756       (((foreach fixnum bignum ratio) (or fixnum bignum ratio))
 757        (multiple-value-bind (q r)
 758            (truncate number divisor)
 759          (values (float q) r)))
 760       (((foreach single-float double-float #!+long-float long-float)
 761         (or rational single-float))
 762        (if (eql divisor 1)
 763            (let ((res (%unary-ftruncate number)))
 764              (values res (- number (coerce res '(dispatch-type number)))))
 765            (ftruncate-float (dispatch-type number))))
 766       #!+long-float
 767       ((long-float (or single-float double-float long-float))
 768        (ftruncate-float long-float))
 769       #!+long-float
 770       (((foreach double-float single-float) long-float)
 771        (ftruncate-float long-float))
 772       ((double-float (or single-float double-float))
 773        (ftruncate-float double-float))
 774       ((single-float double-float)
 775        (ftruncate-float double-float))
 776       (((foreach fixnum bignum ratio)
 777         (foreach single-float double-float #!+long-float long-float))
 778        (ftruncate-float (dispatch-type divisor))))))
 779
 780 (defun ffloor (number &optional (divisor 1))
 781   "Same as FLOOR, but returns first value as a float."
 782   (multiple-value-bind (tru rem) (ftruncate number divisor)
 783     (if (and (not (zerop rem))
 784              (if (minusp divisor)
 785                  (plusp number)
 786                  (minusp number)))
 787         (values (1- tru) (+ rem divisor))
 788         (values tru rem))))
 789
 790 (defun fceiling (number &optional (divisor 1))
 791   "Same as CEILING, but returns first value as a float."
 792   (multiple-value-bind (tru rem) (ftruncate number divisor)
 793     (if (and (not (zerop rem))
 794              (if (minusp divisor)
 795                  (minusp number)
 796                  (plusp number)))
 797         (values (+ tru 1) (- rem divisor))
 798         (values tru rem))))
 799
 800 ;;; FIXME: this probably needs treatment similar to the use of
 801 ;;; %UNARY-FTRUNCATE for FTRUNCATE.
 802 (defun fround (number &optional (divisor 1))
 803   "Same as ROUND, but returns first value as a float."
 804   (multiple-value-bind (res rem)
 805       (round number divisor)
 806     (values (float res (if (floatp rem) rem 1.0)) rem)))
 807 \f
 808 ;;;; comparisons
 809
 810 (defun = (number &rest more-numbers)
 811   #!+sb-doc
 812   "Return T if all of its arguments are numerically equal, NIL otherwise."
 813   (declare (truly-dynamic-extent more-numbers))
 814   (the number number)
 815   (do ((nlist more-numbers (cdr nlist)))
 816       ((atom nlist) t)
 817      (declare (list nlist))
 818      (if (not (= (car nlist) number)) (return nil))))
 819
 820 (defun /= (number &rest more-numbers)
 821   #!+sb-doc
 822   "Return T if no two of its arguments are numerically equal, NIL otherwise."
 823   (declare (truly-dynamic-extent more-numbers))
 824   (do* ((head (the number number) (car nlist))
 825         (nlist more-numbers (cdr nlist)))
 826        ((atom nlist) t)
 827      (declare (list nlist))
 828      (unless (do* ((nl nlist (cdr nl)))
 829                   ((atom nl) t)
 830                (declare (list nl))
 831                (if (= head (car nl)) (return nil)))
 832        (return nil))))
 833
 834 (defun < (number &rest more-numbers)
 835   #!+sb-doc
 836   "Return T if its arguments are in strictly increasing order, NIL otherwise."
 837   (declare (truly-dynamic-extent more-numbers))
 838   (do* ((n (the number number) (car nlist))
 839         (nlist more-numbers (cdr nlist)))
 840        ((atom nlist) t)
 841      (declare (list nlist))
 842      (if (not (< n (car nlist))) (return nil))))
 843
 844 (defun > (number &rest more-numbers)
 845   #!+sb-doc
 846   "Return T if its arguments are in strictly decreasing order, NIL otherwise."
 847   (declare (truly-dynamic-extent more-numbers))
 848   (do* ((n (the number number) (car nlist))
 849         (nlist more-numbers (cdr nlist)))
 850        ((atom nlist) t)
 851      (declare (list nlist))
 852      (if (not (> n (car nlist))) (return nil))))
 853
 854 (defun <= (number &rest more-numbers)
 855   #!+sb-doc
 856   "Return T if arguments are in strictly non-decreasing order, NIL otherwise."
 857   (declare (truly-dynamic-extent more-numbers))
 858   (do* ((n (the number number) (car nlist))
 859         (nlist more-numbers (cdr nlist)))
 860        ((atom nlist) t)
 861      (declare (list nlist))
 862      (if (not (<= n (car nlist))) (return nil))))
 863
 864 (defun >= (number &rest more-numbers)
 865   #!+sb-doc
 866   "Return T if arguments are in strictly non-increasing order, NIL otherwise."
 867   (declare (truly-dynamic-extent more-numbers))
 868   (do* ((n (the number number) (car nlist))
 869         (nlist more-numbers (cdr nlist)))
 870        ((atom nlist) t)
 871      (declare (list nlist))
 872      (if (not (>= n (car nlist))) (return nil))))
 873
 874 (defun max (number &rest more-numbers)
 875   #!+sb-doc
 876   "Return the greatest of its arguments; among EQUALP greatest, return
 877 the first."
 878   (declare (truly-dynamic-extent more-numbers))
 879   (do ((nlist more-numbers (cdr nlist))
 880        (result number))
 881       ((null nlist) (return result))
 882      (declare (list nlist))
 883      (declare (type real number result))
 884      (if (> (car nlist) result) (setq result (car nlist)))))
 885
 886 (defun min (number &rest more-numbers)
 887   #!+sb-doc
 888   "Return the least of its arguments; among EQUALP least, return
 889 the first."
 890   (declare (truly-dynamic-extent more-numbers))
 891   (do ((nlist more-numbers (cdr nlist))
 892        (result number))
 893       ((null nlist) (return result))
 894      (declare (list nlist))
 895      (declare (type real number result))
 896      (if (< (car nlist) result) (setq result (car nlist)))))
 897
 898 (eval-when (:compile-toplevel :execute)
 899
 900 ;;; The INFINITE-X-FINITE-Y and INFINITE-Y-FINITE-X args tell us how
 901 ;;; to handle the case when X or Y is a floating-point infinity and
 902 ;;; the other arg is a rational. (Section 12.1.4.1 of the ANSI spec
 903 ;;; says that comparisons are done by converting the float to a
 904 ;;; rational when comparing with a rational, but infinities can't be
 905 ;;; converted to a rational, so we show some initiative and do it this
 906 ;;; way instead.)
 907 (defun basic-compare (op &key infinite-x-finite-y infinite-y-finite-x)
 908   `(((fixnum fixnum) (,op x y))
 909
 910     ((single-float single-float) (,op x y))
 911     #!+long-float
 912     (((foreach single-float double-float long-float) long-float)
 913      (,op (coerce x 'long-float) y))
 914     #!+long-float
 915     ((long-float (foreach single-float double-float))
 916      (,op x (coerce y 'long-float)))
 917     ((fixnum (foreach single-float double-float))
 918      (if (float-infinity-p y)
 919          ,infinite-y-finite-x
 920          ;; If the fixnum has an exact float representation, do a
 921          ;; float comparison. Otherwise do the slow float -> ratio
 922          ;; conversion.
 923          (multiple-value-bind (lo hi)
 924              (case '(dispatch-type y)
 925                (single-float
 926                 (values most-negative-exactly-single-float-fixnum
 927                         most-positive-exactly-single-float-fixnum))
 928                (double-float
 929                 (values most-negative-exactly-double-float-fixnum
 930                         most-positive-exactly-double-float-fixnum)))
 931            (if (<= lo y hi)
 932                (,op (coerce x '(dispatch-type y)) y)
 933                (,op x (rational y))))))
 934     (((foreach single-float double-float) fixnum)
 935      (if (eql y 0)
 936          (,op x (coerce 0 '(dispatch-type x)))
 937          (if (float-infinity-p x)
 938              ,infinite-x-finite-y
 939              ;; Likewise
 940              (multiple-value-bind (lo hi)
 941                  (case '(dispatch-type x)
 942                    (single-float
 943                     (values most-negative-exactly-single-float-fixnum
 944                             most-positive-exactly-single-float-fixnum))
 945                    (double-float
 946                     (values most-negative-exactly-double-float-fixnum
 947                             most-positive-exactly-double-float-fixnum)))
 948                (if (<= lo y hi)
 949                    (,op x (coerce y '(dispatch-type x)))
 950                    (,op (rational x) y))))))
 951     (((foreach single-float double-float) double-float)
 952      (,op (coerce x 'double-float) y))
 953     ((double-float single-float)
 954      (,op x (coerce y 'double-float)))
 955     (((foreach single-float double-float #!+long-float long-float) rational)
 956      (if (eql y 0)
 957          (,op x (coerce 0 '(dispatch-type x)))
 958          (if (float-infinity-p x)
 959              ,infinite-x-finite-y
 960              (,op (rational x) y))))
 961     (((foreach bignum fixnum ratio) float)
 962      (if (float-infinity-p y)
 963          ,infinite-y-finite-x
 964          (,op x (rational y))))))
 965 ) ; EVAL-WHEN
 966
 967 (macrolet ((def-two-arg-</> (name op ratio-arg1 ratio-arg2 &rest cases)
 968              `(defun ,name (x y)
 969                 (number-dispatch ((x real) (y real))
 970                                  (basic-compare
 971                                   ,op
 972                                   :infinite-x-finite-y
 973                                   (,op x (coerce 0 '(dispatch-type x)))
 974                                   :infinite-y-finite-x
 975                                   (,op (coerce 0 '(dispatch-type y)) y))
 976                                  (((foreach fixnum bignum) ratio)
 977                                   (,op x (,ratio-arg2 (numerator y)
 978                                                       (denominator y))))
 979                                  ((ratio integer)
 980                                   (,op (,ratio-arg1 (numerator x)
 981                                                     (denominator x))
 982                                        y))
 983                                  ((ratio ratio)
 984                                   (,op (* (numerator   (truly-the ratio x))
 985                                           (denominator (truly-the ratio y)))
 986                                        (* (numerator   (truly-the ratio y))
 987                                           (denominator (truly-the ratio x)))))
 988                                  ,@cases))))
 989   (def-two-arg-</> two-arg-< < floor ceiling
 990     ((fixnum bignum)
 991      (bignum-plus-p y))
 992     ((bignum fixnum)
 993      (not (bignum-plus-p x)))
 994     ((bignum bignum)
 995      (minusp (bignum-compare x y))))
 996   (def-two-arg-</> two-arg-> > ceiling floor
 997     ((fixnum bignum)
 998      (not (bignum-plus-p y)))
 999     ((bignum fixnum)
1000      (bignum-plus-p x))
1001     ((bignum bignum)
1002      (plusp (bignum-compare x y)))))
1003
1004 (defun two-arg-= (x y)
1005   (number-dispatch ((x number) (y number))
1006     (basic-compare =
1007                    ;; An infinite value is never equal to a finite value.
1008                    :infinite-x-finite-y nil
1009                    :infinite-y-finite-x nil)
1010     ((fixnum (or bignum ratio)) nil)
1011
1012     ((bignum (or fixnum ratio)) nil)
1013     ((bignum bignum)
1014      (zerop (bignum-compare x y)))
1015
1016     ((ratio integer) nil)
1017     ((ratio ratio)
1018      (and (eql (numerator x) (numerator y))
1019           (eql (denominator x) (denominator y))))
1020
1021     ((complex complex)
1022      (and (= (realpart x) (realpart y))
1023           (= (imagpart x) (imagpart y))))
1024     (((foreach fixnum bignum ratio single-float double-float
1025                #!+long-float long-float) complex)
1026      (and (= x (realpart y))
1027           (zerop (imagpart y))))
1028     ((complex (or float rational))
1029      (and (= (realpart x) y)
1030           (zerop (imagpart x))))))
1031 \f
1032 ;;;; logicals
1033
1034 (defun logior (&rest integers)
1035   #!+sb-doc
1036   "Return the bit-wise or of its arguments. Args must be integers."
1037   (declare (list integers))
1038   (declare (truly-dynamic-extent integers))
1039   (if integers
1040       (do ((result (pop integers) (logior result (pop integers))))
1041           ((null integers) result)
1042         (declare (integer result)))
1043       0))
1044
1045 (defun logxor (&rest integers)
1046   #!+sb-doc
1047   "Return the bit-wise exclusive or of its arguments. Args must be integers."
1048   (declare (list integers))
1049   (declare (truly-dynamic-extent integers))
1050   (if integers
1051       (do ((result (pop integers) (logxor result (pop integers))))
1052           ((null integers) result)
1053         (declare (integer result)))
1054       0))
1055
1056 (defun logand (&rest integers)
1057   #!+sb-doc
1058   "Return the bit-wise and of its arguments. Args must be integers."
1059   (declare (list integers))
1060   (declare (truly-dynamic-extent integers))
1061   (if integers
1062       (do ((result (pop integers) (logand result (pop integers))))
1063           ((null integers) result)
1064         (declare (integer result)))
1065       -1))
1066
1067 (defun logeqv (&rest integers)
1068   #!+sb-doc
1069   "Return the bit-wise equivalence of its arguments. Args must be integers."
1070   (declare (list integers))
1071   (declare (truly-dynamic-extent integers))
1072   (if integers
1073       (do ((result (pop integers) (logeqv result (pop integers))))
1074           ((null integers) result)
1075         (declare (integer result)))
1076       -1))
1077
1078 (defun lognot (number)
1079   #!+sb-doc
1080   "Return the bit-wise logical not of integer."
1081   (etypecase number
1082     (fixnum (lognot (truly-the fixnum number)))
1083     (bignum (bignum-logical-not number))))
1084
1085 (macrolet ((def (name op big-op &optional doc)
1086              `(defun ,name (integer1 integer2)
1087                 ,@(when doc
1088                     (list doc))
1089                 (let ((x integer1)
1090                       (y integer2))
1091                   (number-dispatch ((x integer) (y integer))
1092                     (bignum-cross-fixnum ,op ,big-op))))))
1093   (def two-arg-and logand bignum-logical-and)
1094   (def two-arg-ior logior bignum-logical-ior)
1095   (def two-arg-xor logxor bignum-logical-xor)
1096   ;; BIGNUM-LOGICAL-{AND,IOR,XOR} need not return a bignum, so must
1097   ;; call the generic LOGNOT...
1098   (def two-arg-eqv logeqv (lambda (x y) (lognot (bignum-logical-xor x y))))
1099   (def lognand lognand
1100        (lambda (x y) (lognot (bignum-logical-and x y)))
1101        #!+sb-doc "Complement the logical AND of INTEGER1 and INTEGER2.")
1102   (def lognor lognor
1103        (lambda (x y) (lognot (bignum-logical-ior x y)))
1104        #!+sb-doc "Complement the logical AND of INTEGER1 and INTEGER2.")
1105   ;; ... but BIGNUM-LOGICAL-NOT on a bignum will always return a bignum
1106   (def logandc1 logandc1
1107        (lambda (x y) (bignum-logical-and (bignum-logical-not x) y))
1108        #!+sb-doc "Bitwise AND (LOGNOT INTEGER1) with INTEGER2.")
1109   (def logandc2 logandc2
1110        (lambda (x y) (bignum-logical-and x (bignum-logical-not y)))
1111        #!+sb-doc "Bitwise AND INTEGER1 with (LOGNOT INTEGER2).")
1112   (def logorc1 logorc1
1113        (lambda (x y) (bignum-logical-ior (bignum-logical-not x) y))
1114        #!+sb-doc "Bitwise OR (LOGNOT INTEGER1) with INTEGER2.")
1115   (def logorc2 logorc2
1116        (lambda (x y) (bignum-logical-ior x (bignum-logical-not y)))
1117        #!+sb-doc "Bitwise OR INTEGER1 with (LOGNOT INTEGER2)."))
1118
1119 (defun logcount (integer)
1120   #!+sb-doc
1121   "Count the number of 1 bits if INTEGER is positive, and the number of 0 bits
1122   if INTEGER is negative."
1123   (etypecase integer
1124     (fixnum
1125      (logcount (truly-the (integer 0
1126                                    #.(max sb!xc:most-positive-fixnum
1127                                           (lognot sb!xc:most-negative-fixnum)))
1128                           (if (minusp (truly-the fixnum integer))
1129                               (lognot (truly-the fixnum integer))
1130                               integer))))
1131     (bignum
1132      (bignum-logcount integer))))
1133
1134 (defun logtest (integer1 integer2)
1135   #!+sb-doc
1136   "Predicate which returns T if logand of integer1 and integer2 is not zero."
1137   (logtest integer1 integer2))
1138
1139 (defun logbitp (index integer)
1140   #!+sb-doc
1141   "Predicate returns T if bit index of integer is a 1."
1142   (number-dispatch ((index integer) (integer integer))
1143     ((fixnum fixnum) (if (< index sb!vm:n-positive-fixnum-bits)
1144                          (not (zerop (logand integer (ash 1 index))))
1145                          (minusp integer)))
1146     ((fixnum bignum) (bignum-logbitp index integer))
1147     ((bignum (foreach fixnum bignum)) (minusp integer))))
1148
1149 (defun ash (integer count)
1150   #!+sb-doc
1151   "Shifts integer left by count places preserving sign. - count shifts right."
1152   (declare (integer integer count))
1153   (etypecase integer
1154     (fixnum
1155      (cond ((zerop integer)
1156             0)
1157            ((fixnump count)
1158             (let ((length (integer-length (truly-the fixnum integer)))
1159                   (count (truly-the fixnum count)))
1160               (declare (fixnum length count))
1161               (cond ((and (plusp count)
1162                           (> (+ length count)
1163                              (integer-length most-positive-fixnum)))
1164                      (bignum-ashift-left (make-small-bignum integer) count))
1165                     (t
1166                      (truly-the fixnum
1167                                 (ash (truly-the fixnum integer) count))))))
1168            ((minusp count)
1169             (if (minusp integer) -1 0))
1170            (t
1171             (bignum-ashift-left (make-small-bignum integer) count))))
1172     (bignum
1173      (if (plusp count)
1174          (bignum-ashift-left integer count)
1175          (bignum-ashift-right integer (- count))))))
1176
1177 (defun integer-length (integer)
1178   #!+sb-doc
1179   "Return the number of non-sign bits in the twos-complement representation
1180   of INTEGER."
1181   (etypecase integer
1182     (fixnum
1183      (integer-length (truly-the fixnum integer)))
1184     (bignum
1185      (bignum-integer-length integer))))
1186 \f
1187 ;;;; BYTE, bytespecs, and related operations
1188
1189 (defun byte (size position)
1190   #!+sb-doc
1191   "Return a byte specifier which may be used by other byte functions
1192   (e.g. LDB)."
1193   (byte size position))
1194
1195 (defun byte-size (bytespec)
1196   #!+sb-doc
1197   "Return the size part of the byte specifier bytespec."
1198   (byte-size bytespec))
1199
1200 (defun byte-position (bytespec)
1201   #!+sb-doc
1202   "Return the position part of the byte specifier bytespec."
1203   (byte-position bytespec))
1204
1205 (defun ldb (bytespec integer)
1206   #!+sb-doc
1207   "Extract the specified byte from integer, and right justify result."
1208   (ldb bytespec integer))
1209
1210 (defun ldb-test (bytespec integer)
1211   #!+sb-doc
1212   "Return T if any of the specified bits in integer are 1's."
1213   (ldb-test bytespec integer))
1214
1215 (defun mask-field (bytespec integer)
1216   #!+sb-doc
1217   "Extract the specified byte from integer,  but do not right justify result."
1218   (mask-field bytespec integer))
1219
1220 (defun dpb (newbyte bytespec integer)
1221   #!+sb-doc
1222   "Return new integer with newbyte in specified position, newbyte is right justified."
1223   (dpb newbyte bytespec integer))
1224
1225 (defun deposit-field (newbyte bytespec integer)
1226   #!+sb-doc
1227   "Return new integer with newbyte in specified position, newbyte is not right justified."
1228   (deposit-field newbyte bytespec integer))
1229
1230 (defun %ldb (size posn integer)
1231   (logand (ash integer (- posn))
1232           (1- (ash 1 size))))
1233
1234 (defun %mask-field (size posn integer)
1235   (logand integer (ash (1- (ash 1 size)) posn)))
1236
1237 (defun %dpb (newbyte size posn integer)
1238   (let ((mask (1- (ash 1 size))))
1239     (logior (logand integer (lognot (ash mask posn)))
1240             (ash (logand newbyte mask) posn))))
1241
1242 (defun %deposit-field (newbyte size posn integer)
1243   (let ((mask (ash (ldb (byte size 0) -1) posn)))
1244     (logior (logand newbyte mask)
1245             (logand integer (lognot mask)))))
1246
1247 (defun sb!c::mask-signed-field (size integer)
1248   #!+sb-doc
1249   "Extract SIZE lower bits from INTEGER, considering them as a
1250 2-complement SIZE-bits representation of a signed integer."
1251   (cond ((zerop size)
1252          0)
1253         ((logbitp (1- size) integer)
1254          (dpb integer (byte size 0) -1))
1255         (t
1256          (ldb (byte size 0) integer))))
1257
1258 \f
1259 ;;;; BOOLE
1260
1261 ;;; The boole function dispaches to any logic operation depending on
1262 ;;;     the value of a variable. Presently, legal selector values are [0..15].
1263 ;;;     boole is open coded for calls with a constant selector. or with calls
1264 ;;;     using any of the constants declared below.
1265
1266 (defconstant boole-clr 0
1267   #!+sb-doc
1268   "Boole function op, makes BOOLE return 0.")
1269
1270 (defconstant boole-set 1
1271   #!+sb-doc
1272   "Boole function op, makes BOOLE return -1.")
1273
1274 (defconstant boole-1   2
1275   #!+sb-doc
1276   "Boole function op, makes BOOLE return integer1.")
1277
1278 (defconstant boole-2   3
1279   #!+sb-doc
1280   "Boole function op, makes BOOLE return integer2.")
1281
1282 (defconstant boole-c1  4
1283   #!+sb-doc
1284   "Boole function op, makes BOOLE return complement of integer1.")
1285
1286 (defconstant boole-c2  5
1287   #!+sb-doc
1288   "Boole function op, makes BOOLE return complement of integer2.")
1289
1290 (defconstant boole-and 6
1291   #!+sb-doc
1292   "Boole function op, makes BOOLE return logand of integer1 and integer2.")
1293
1294 (defconstant boole-ior 7
1295   #!+sb-doc
1296   "Boole function op, makes BOOLE return logior of integer1 and integer2.")
1297
1298 (defconstant boole-xor 8
1299   #!+sb-doc
1300   "Boole function op, makes BOOLE return logxor of integer1 and integer2.")
1301
1302 (defconstant boole-eqv 9
1303   #!+sb-doc
1304   "Boole function op, makes BOOLE return logeqv of integer1 and integer2.")
1305
1306 (defconstant boole-nand  10
1307   #!+sb-doc
1308   "Boole function op, makes BOOLE return log nand of integer1 and integer2.")
1309
1310 (defconstant boole-nor   11
1311   #!+sb-doc
1312   "Boole function op, makes BOOLE return lognor of integer1 and integer2.")
1313
1314 (defconstant boole-andc1 12
1315   #!+sb-doc
1316   "Boole function op, makes BOOLE return logandc1 of integer1 and integer2.")
1317
1318 (defconstant boole-andc2 13
1319   #!+sb-doc
1320   "Boole function op, makes BOOLE return logandc2 of integer1 and integer2.")
1321
1322 (defconstant boole-orc1  14
1323   #!+sb-doc
1324   "Boole function op, makes BOOLE return logorc1 of integer1 and integer2.")
1325
1326 (defconstant boole-orc2  15
1327   #!+sb-doc
1328   "Boole function op, makes BOOLE return logorc2 of integer1 and integer2.")
1329
1330 (defun boole (op integer1 integer2)
1331   #!+sb-doc
1332   "Bit-wise boolean function on two integers. Function chosen by OP:
1333         0       BOOLE-CLR
1334         1       BOOLE-SET
1335         2       BOOLE-1
1336         3       BOOLE-2
1337         4       BOOLE-C1
1338         5       BOOLE-C2
1339         6       BOOLE-AND
1340         7       BOOLE-IOR
1341         8       BOOLE-XOR
1342         9       BOOLE-EQV
1343         10      BOOLE-NAND
1344         11      BOOLE-NOR
1345         12      BOOLE-ANDC1
1346         13      BOOLE-ANDC2
1347         14      BOOLE-ORC1
1348         15      BOOLE-ORC2"
1349   (case op
1350     (0 (boole 0 integer1 integer2))
1351     (1 (boole 1 integer1 integer2))
1352     (2 (boole 2 integer1 integer2))
1353     (3 (boole 3 integer1 integer2))
1354     (4 (boole 4 integer1 integer2))
1355     (5 (boole 5 integer1 integer2))
1356     (6 (boole 6 integer1 integer2))
1357     (7 (boole 7 integer1 integer2))
1358     (8 (boole 8 integer1 integer2))
1359     (9 (boole 9 integer1 integer2))
1360     (10 (boole 10 integer1 integer2))
1361     (11 (boole 11 integer1 integer2))
1362     (12 (boole 12 integer1 integer2))
1363     (13 (boole 13 integer1 integer2))
1364     (14 (boole 14 integer1 integer2))
1365     (15 (boole 15 integer1 integer2))
1366     (t (error 'type-error :datum op :expected-type '(mod 16)))))
1367 \f
1368 ;;;; GCD and LCM
1369
1370 (defun gcd (&rest integers)
1371   #!+sb-doc
1372   "Return the greatest common divisor of the arguments, which must be
1373   integers. Gcd with no arguments is defined to be 0."
1374   (declare (truly-dynamic-extent integers))
1375   (cond ((null integers) 0)
1376         ((null (cdr integers)) (abs (the integer (car integers))))
1377         (t
1378          (do ((gcd (the integer (car integers))
1379                    (gcd gcd (the integer (car rest))))
1380               (rest (cdr integers) (cdr rest)))
1381              ((null rest) gcd)
1382            (declare (integer gcd)
1383                     (list rest))))))
1384
1385 (defun lcm (&rest integers)
1386   #!+sb-doc
1387   "Return the least common multiple of one or more integers. LCM of no
1388   arguments is defined to be 1."
1389   (declare (truly-dynamic-extent integers))
1390   (cond ((null integers) 1)
1391         ((null (cdr integers)) (abs (the integer (car integers))))
1392         (t
1393          (do ((lcm (the integer (car integers))
1394                    (lcm lcm (the integer (car rest))))
1395               (rest (cdr integers) (cdr rest)))
1396              ((null rest) lcm)
1397            (declare (integer lcm) (list rest))))))
1398
1399 (defun two-arg-lcm (n m)
1400   (declare (integer n m))
1401   (if (or (zerop n) (zerop m))
1402       0
1403       ;; KLUDGE: I'm going to assume that it was written this way
1404       ;; originally for a reason.  However, this is a somewhat
1405       ;; complicated way of writing the algorithm in the CLHS page for
1406       ;; LCM, and I don't know why.  To be investigated.  -- CSR,
1407       ;; 2003-09-11
1408       ;;
1409       ;;    It seems to me that this is written this way to avoid
1410       ;;    unnecessary bignumification of intermediate results.
1411       ;;        -- TCR, 2008-03-05
1412       (let ((m (abs m))
1413             (n (abs n)))
1414         (multiple-value-bind (max min)
1415             (if (> m n)
1416                 (values m n)
1417                 (values n m))
1418           (* (truncate max (gcd n m)) min)))))
1419
1420 ;;; Do the GCD of two integer arguments. With fixnum arguments, we use the
1421 ;;; binary GCD algorithm from Knuth's seminumerical algorithms (slightly
1422 ;;; structurified), otherwise we call BIGNUM-GCD. We pick off the special case
1423 ;;; of 0 before the dispatch so that the bignum code doesn't have to worry
1424 ;;; about "small bignum" zeros.
1425 (defun two-arg-gcd (u v)
1426   (cond ((eql u 0) (abs v))
1427         ((eql v 0) (abs u))
1428         (t
1429          (number-dispatch ((u integer) (v integer))
1430            ((fixnum fixnum)
1431             (locally
1432                 (declare (optimize (speed 3) (safety 0)))
1433               (do ((k 0 (1+ k))
1434                    (u (abs u) (ash u -1))
1435                    (v (abs v) (ash v -1)))
1436                   ((oddp (logior u v))
1437                    (do ((temp (if (oddp u) (- v) (ash u -1))
1438                               (ash temp -1)))
1439                        (nil)
1440                      (declare (fixnum temp))
1441                      (when (oddp temp)
1442                        (if (plusp temp)
1443                            (setq u temp)
1444                            (setq v (- temp)))
1445                        (setq temp (- u v))
1446                        (when (zerop temp)
1447                          (let ((res (ash u k)))
1448                            (declare (type sb!vm:signed-word res)
1449                                     (optimize (inhibit-warnings 3)))
1450                            (return res))))))
1451                 (declare (type (mod #.sb!vm:n-word-bits) k)
1452                          (type sb!vm:signed-word u v)))))
1453            ((bignum bignum)
1454             (bignum-gcd u v))
1455            ((bignum fixnum)
1456             (bignum-gcd u (make-small-bignum v)))
1457            ((fixnum bignum)
1458             (bignum-gcd (make-small-bignum u) v))))))
1459 \f
1460 ;;;; from Robert Smith; slightly changed not to cons unnecessarily.
1461 (defun isqrt (n)
1462   #!+sb-doc
1463   "Return the greatest integer less than or equal to the square root of N."
1464   (declare (type unsigned-byte n))
1465   (cond
1466     ((> n 24)
1467      (let* ((n-fourth-size (ash (1- (integer-length n)) -2))
1468             (n-significant-half (ash n (- (ash n-fourth-size 1))))
1469             (n-significant-half-isqrt (isqrt n-significant-half))
1470             (zeroth-iteration (ash n-significant-half-isqrt n-fourth-size)))
1471        (multiple-value-bind (quot rem)
1472            (floor n zeroth-iteration)
1473          (let ((first-iteration (ash (+ zeroth-iteration quot) -1)))
1474            (cond ((oddp quot)
1475                   first-iteration)
1476                  ((> (expt (- first-iteration zeroth-iteration) 2) rem)
1477                   (1- first-iteration))
1478                  (t
1479                   first-iteration))))))
1480     ((> n 15) 4)
1481     ((> n  8) 3)
1482     ((> n  3) 2)
1483     ((> n  0) 1)
1484     ((= n  0) 0)))
1485 \f
1486 ;;;; miscellaneous number predicates
1487
1488 (macrolet ((def (name doc)
1489              `(defun ,name (number) ,doc (,name number))))
1490   (def zerop "Is this number zero?")
1491   (def plusp "Is this real number strictly positive?")
1492   (def minusp "Is this real number strictly negative?")
1493   (def oddp "Is this integer odd?")
1494   (def evenp "Is this integer even?"))
1495 \f
1496 ;;;; modular functions
1497 #.
1498 (collect ((forms))
1499   (flet ((unsigned-definition (name lambda-list width)
1500            (let ((pattern (1- (ash 1 width))))
1501              `(defun ,name ,lambda-list
1502                (flet ((prepare-argument (x)
1503                         (declare (integer x))
1504                         (etypecase x
1505                           ((unsigned-byte ,width) x)
1506                           (fixnum (logand x ,pattern))
1507                           (bignum (logand x ,pattern)))))
1508                  (,name ,@(loop for arg in lambda-list
1509                                 collect `(prepare-argument ,arg)))))))
1510          (signed-definition (name lambda-list width)
1511            `(defun ,name ,lambda-list
1512               (flet ((prepare-argument (x)
1513                        (declare (integer x))
1514                        (etypecase x
1515                          ((signed-byte ,width) x)
1516                          (fixnum (sb!c::mask-signed-field ,width x))
1517                          (bignum (sb!c::mask-signed-field ,width x)))))
1518                 (,name ,@(loop for arg in lambda-list
1519                                collect `(prepare-argument ,arg)))))))
1520     (flet ((do-mfuns (class)
1521              (loop for infos being each hash-value of (sb!c::modular-class-funs class)
1522                    ;; FIXME: We need to process only "toplevel" functions
1523                    when (listp infos)
1524                    do (loop for info in infos
1525                             for name = (sb!c::modular-fun-info-name info)
1526                             and width = (sb!c::modular-fun-info-width info)
1527                             and signedp = (sb!c::modular-fun-info-signedp info)
1528                             and lambda-list = (sb!c::modular-fun-info-lambda-list info)
1529                             if signedp
1530                             do (forms (signed-definition name lambda-list width))
1531                             else
1532                             do (forms (unsigned-definition name lambda-list width))))))
1533       (do-mfuns sb!c::*untagged-unsigned-modular-class*)
1534       (do-mfuns sb!c::*untagged-signed-modular-class*)
1535       (do-mfuns sb!c::*tagged-modular-class*)))
1536   `(progn ,@(sort (forms) #'string< :key #'cadr)))
1537
1538 ;;; KLUDGE: these out-of-line definitions can't use the modular
1539 ;;; arithmetic, as that is only (currently) defined for constant
1540 ;;; shifts.  See also the comment in (LOGAND OPTIMIZER) for more
1541 ;;; discussion of this hack.  -- CSR, 2003-10-09
1542 #!+#.(cl:if (cl:= sb!vm:n-machine-word-bits 32) '(and) '(or))
1543 (defun sb!vm::ash-left-mod32 (integer amount)
1544   (etypecase integer
1545     ((unsigned-byte 32) (ldb (byte 32 0) (ash integer amount)))
1546     (fixnum (ldb (byte 32 0) (ash (logand integer #xffffffff) amount)))
1547     (bignum (ldb (byte 32 0) (ash (logand integer #xffffffff) amount)))))
1548 #!+#.(cl:if (cl:= sb!vm:n-machine-word-bits 64) '(and) '(or))
1549 (defun sb!vm::ash-left-mod64 (integer amount)
1550   (etypecase integer
1551     ((unsigned-byte 64) (ldb (byte 64 0) (ash integer amount)))
1552     (fixnum (ldb (byte 64 0) (ash (logand integer #xffffffffffffffff) amount)))
1553     (bignum (ldb (byte 64 0)
1554                  (ash (logand integer #xffffffffffffffff) amount)))))
1555
1556 #!+(or x86 x86-64)
1557 (defun sb!vm::ash-left-modfx (integer amount)
1558   (let ((fixnum-width (- sb!vm:n-word-bits sb!vm:n-fixnum-tag-bits)))
1559     (etypecase integer
1560       (fixnum (sb!c::mask-signed-field fixnum-width (ash integer amount)))
1561       (integer (sb!c::mask-signed-field fixnum-width (ash (sb!c::mask-signed-field fixnum-width integer) amount))))))