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