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 upgraded-complex-part-type (spec)
 215   #!+sb-doc
 216   "Return the element type of the most specialized COMPLEX number type that
 217    can hold parts of type SPEC."
 218   (cond ((unknown-type-p (specifier-type spec))
 219          (error "undefined type: ~S" spec))
 220         ((subtypep spec 'single-float)
 221          'single-float)
 222         ((subtypep spec 'double-float)
 223          'double-float)
 224         #!+long-float
 225         ((subtypep spec 'long-float)
 226          'long-float)
 227         ((subtypep spec 'rational)
 228          'rational)
 229         (t
 230          'real)))
 231
 232 (defun complex (realpart &optional (imagpart 0))
 233   #!+sb-doc
 234   "Return a complex number with the specified real and imaginary components."
 235   (flet ((%%make-complex (realpart imagpart)
 236            (cond #!+long-float
 237                  ((and (typep realpart 'long-float)
 238                        (typep imagpart 'long-float))
 239                   (truly-the (complex long-float)
 240                              (complex realpart imagpart)))
 241                  ((and (typep realpart 'double-float)
 242                        (typep imagpart 'double-float))
 243                   (truly-the (complex double-float)
 244                              (complex realpart imagpart)))
 245                  ((and (typep realpart 'single-float)
 246                        (typep imagpart 'single-float))
 247                   (truly-the (complex single-float)
 248                              (complex realpart imagpart)))
 249                  (t
 250                   (%make-complex realpart imagpart)))))
 251   (number-dispatch ((realpart real) (imagpart real))
 252     ((rational rational)
 253      (canonical-complex realpart imagpart))
 254     (float-contagion %%make-complex realpart imagpart (rational)))))
 255
 256 (defun realpart (number)
 257   #!+sb-doc
 258   "Extract the real part of a number."
 259   (typecase number
 260     #!+long-float
 261     ((complex long-float)
 262      (truly-the long-float (realpart number)))
 263     ((complex double-float)
 264      (truly-the double-float (realpart number)))
 265     ((complex single-float)
 266      (truly-the single-float (realpart number)))
 267     ((complex rational)
 268      (sb!kernel:%realpart number))
 269     (t
 270      number)))
 271
 272 (defun imagpart (number)
 273   #!+sb-doc
 274   "Extract the imaginary part of a number."
 275   (typecase number
 276     #!+long-float
 277     ((complex long-float)
 278      (truly-the long-float (imagpart number)))
 279     ((complex double-float)
 280      (truly-the double-float (imagpart number)))
 281     ((complex single-float)
 282      (truly-the single-float (imagpart number)))
 283     ((complex rational)
 284      (sb!kernel:%imagpart number))
 285     (float
 286      (float 0 number))
 287     (t
 288      0)))
 289
 290 (defun conjugate (number)
 291   #!+sb-doc
 292   "Return the complex conjugate of NUMBER. For non-complex numbers, this is
 293   an identity."
 294   (if (complexp number)
 295       (complex (realpart number) (- (imagpart number)))
 296       number))
 297
 298 (defun signum (number)
 299   #!+sb-doc
 300   "If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))."
 301   (if (zerop number)
 302       number
 303       (if (rationalp number)
 304           (if (plusp number) 1 -1)
 305           (/ number (abs number)))))
 306 \f
 307 ;;;; ratios
 308
 309 (defun numerator (number)
 310   #!+sb-doc
 311   "Return the numerator of NUMBER, which must be rational."
 312   (numerator number))
 313
 314 (defun denominator (number)
 315   #!+sb-doc
 316   "Return the denominator of NUMBER, which must be rational."
 317   (denominator number))
 318 \f
 319 ;;;; arithmetic operations
 320
 321 (macrolet ((define-arith (op init doc)
 322              #!-sb-doc (declare (ignore doc))
 323              `(defun ,op (&rest args)
 324                 #!+sb-doc ,doc
 325                 (if (null args) ,init
 326                     (do ((args (cdr args) (cdr args))
 327                          (result (car args) (,op result (car args))))
 328                         ((null args) result)
 329                       ;; to signal TYPE-ERROR when exactly 1 arg of wrong type:
 330                       (declare (type number result)))))))
 331   (define-arith + 0
 332     "Return the sum of its arguments. With no args, returns 0.")
 333   (define-arith * 1
 334     "Return the product of its arguments. With no args, returns 1."))
 335
 336 (defun - (number &rest more-numbers)
 337   #!+sb-doc
 338   "Subtract the second and all subsequent arguments from the first;
 339   or with one argument, negate the first argument."
 340   (if more-numbers
 341       (do ((nlist more-numbers (cdr nlist))
 342            (result number))
 343           ((atom nlist) result)
 344          (declare (list nlist))
 345          (setq result (- result (car nlist))))
 346       (- number)))
 347
 348 (defun / (number &rest more-numbers)
 349   #!+sb-doc
 350   "Divide the first argument by each of the following arguments, in turn.
 351   With one argument, return reciprocal."
 352   (if more-numbers
 353       (do ((nlist more-numbers (cdr nlist))
 354            (result number))
 355           ((atom nlist) result)
 356          (declare (list nlist))
 357          (setq result (/ result (car nlist))))
 358       (/ number)))
 359
 360 (defun 1+ (number)
 361   #!+sb-doc
 362   "Return NUMBER + 1."
 363   (1+ number))
 364
 365 (defun 1- (number)
 366   #!+sb-doc
 367   "Return NUMBER - 1."
 368   (1- number))
 369
 370 (eval-when (:compile-toplevel)
 371
 372 (sb!xc:defmacro two-arg-+/- (name op big-op)
 373   `(defun ,name (x y)
 374      (number-dispatch ((x number) (y number))
 375        (bignum-cross-fixnum ,op ,big-op)
 376        (float-contagion ,op x y)
 377
 378        ((complex complex)
 379         (canonical-complex (,op (realpart x) (realpart y))
 380                            (,op (imagpart x) (imagpart y))))
 381        (((foreach bignum fixnum ratio single-float double-float
 382                   #!+long-float long-float) complex)
 383         (complex (,op x (realpart y)) (,op (imagpart y))))
 384        ((complex (or rational float))
 385         (complex (,op (realpart x) y) (imagpart x)))
 386
 387        (((foreach fixnum bignum) ratio)
 388         (let* ((dy (denominator y))
 389                (n (,op (* x dy) (numerator y))))
 390           (%make-ratio n dy)))
 391        ((ratio integer)
 392         (let* ((dx (denominator x))
 393                (n (,op (numerator x) (* y dx))))
 394           (%make-ratio n dx)))
 395        ((ratio ratio)
 396         (let* ((nx (numerator x))
 397                (dx (denominator x))
 398                (ny (numerator y))
 399                (dy (denominator y))
 400                (g1 (gcd dx dy)))
 401           (if (eql g1 1)
 402               (%make-ratio (,op (* nx dy) (* dx ny)) (* dx dy))
 403               (let* ((t1 (,op (* nx (truncate dy g1)) (* (truncate dx g1) ny)))
 404                      (g2 (gcd t1 g1))
 405                      (t2 (truncate dx g1)))
 406                 (cond ((eql t1 0) 0)
 407                       ((eql g2 1)
 408                        (%make-ratio t1 (* t2 dy)))
 409                       (T (let* ((nn (truncate t1 g2))
 410                                 (t3 (truncate dy g2))
 411                                 (nd (if (eql t2 1) t3 (* t2 t3))))
 412                            (if (eql nd 1) nn (%make-ratio nn nd))))))))))))
 413
 414 ) ; EVAL-WHEN
 415
 416 (two-arg-+/- two-arg-+ + add-bignums)
 417 (two-arg-+/- two-arg-- - subtract-bignum)
 418
 419 (defun two-arg-* (x y)
 420   (flet ((integer*ratio (x y)
 421            (if (eql x 0) 0
 422                (let* ((ny (numerator y))
 423                       (dy (denominator y))
 424                       (gcd (gcd x dy)))
 425                  (if (eql gcd 1)
 426                      (%make-ratio (* x ny) dy)
 427                      (let ((nn (* (truncate x gcd) ny))
 428                            (nd (truncate dy gcd)))
 429                        (if (eql nd 1)
 430                            nn
 431                            (%make-ratio nn nd)))))))
 432          (complex*real (x y)
 433            (canonical-complex (* (realpart x) y) (* (imagpart x) y))))
 434     (number-dispatch ((x number) (y number))
 435       (float-contagion * x y)
 436
 437       ((fixnum fixnum) (multiply-fixnums x y))
 438       ((bignum fixnum) (multiply-bignum-and-fixnum x y))
 439       ((fixnum bignum) (multiply-bignum-and-fixnum y x))
 440       ((bignum bignum) (multiply-bignums x y))
 441
 442       ((complex complex)
 443        (let* ((rx (realpart x))
 444               (ix (imagpart x))
 445               (ry (realpart y))
 446               (iy (imagpart y)))
 447          (canonical-complex (- (* rx ry) (* ix iy)) (+ (* rx iy) (* ix ry)))))
 448       (((foreach bignum fixnum ratio single-float double-float
 449                  #!+long-float long-float)
 450         complex)
 451        (complex*real y x))
 452       ((complex (or rational float))
 453        (complex*real x y))
 454
 455       (((foreach bignum fixnum) ratio) (integer*ratio x y))
 456       ((ratio integer) (integer*ratio y x))
 457       ((ratio ratio)
 458        (let* ((nx (numerator x))
 459               (dx (denominator x))
 460               (ny (numerator y))
 461               (dy (denominator y))
 462               (g1 (gcd nx dy))
 463               (g2 (gcd dx ny)))
 464          (build-ratio (* (maybe-truncate nx g1)
 465                          (maybe-truncate ny g2))
 466                       (* (maybe-truncate dx g2)
 467                          (maybe-truncate dy g1))))))))
 468
 469 ;;; Divide two integers, producing a canonical rational. If a fixnum,
 470 ;;; we see whether they divide evenly before trying the GCD. In the
 471 ;;; bignum case, we don't bother, since bignum division is expensive,
 472 ;;; and the test is not very likely to succeed.
 473 (defun integer-/-integer (x y)
 474   (if (and (typep x 'fixnum) (typep y 'fixnum))
 475       (multiple-value-bind (quo rem) (truncate x y)
 476         (if (zerop rem)
 477             quo
 478             (let ((gcd (gcd x y)))
 479               (declare (fixnum gcd))
 480               (if (eql gcd 1)
 481                   (build-ratio x y)
 482                   (build-ratio (truncate x gcd) (truncate y gcd))))))
 483       (let ((gcd (gcd x y)))
 484         (if (eql gcd 1)
 485             (build-ratio x y)
 486             (build-ratio (truncate x gcd) (truncate y gcd))))))
 487
 488 (defun two-arg-/ (x y)
 489   (number-dispatch ((x number) (y number))
 490     (float-contagion / x y (ratio integer))
 491
 492     ((complex complex)
 493      (let* ((rx (realpart x))
 494             (ix (imagpart x))
 495             (ry (realpart y))
 496             (iy (imagpart y)))
 497        (if (> (abs ry) (abs iy))
 498            (let* ((r (/ iy ry))
 499                   (dn (* ry (+ 1 (* r r)))))
 500              (canonical-complex (/ (+ rx (* ix r)) dn)
 501                                 (/ (- ix (* rx r)) dn)))
 502            (let* ((r (/ ry iy))
 503                   (dn (* iy (+ 1 (* r r)))))
 504              (canonical-complex (/ (+ (* rx r) ix) dn)
 505                                 (/ (- (* ix r) rx) dn))))))
 506     (((foreach integer ratio single-float double-float) complex)
 507      (let* ((ry (realpart y))
 508             (iy (imagpart y)))
 509        (if (> (abs ry) (abs iy))
 510            (let* ((r (/ iy ry))
 511                   (dn (* ry (+ 1 (* r r)))))
 512              (canonical-complex (/ x dn)
 513                                 (/ (- (* x r)) dn)))
 514            (let* ((r (/ ry iy))
 515                   (dn (* iy (+ 1 (* r r)))))
 516              (canonical-complex (/ (* x r) dn)
 517                                 (/ (- x) dn))))))
 518     ((complex (or rational float))
 519      (canonical-complex (/ (realpart x) y)
 520                         (/ (imagpart x) y)))
 521
 522     ((ratio ratio)
 523      (let* ((nx (numerator x))
 524             (dx (denominator x))
 525             (ny (numerator y))
 526             (dy (denominator y))
 527             (g1 (gcd nx ny))
 528             (g2 (gcd dx dy)))
 529        (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2))
 530                     (* (maybe-truncate dx g2) (maybe-truncate ny g1)))))
 531
 532     ((integer integer)
 533      (integer-/-integer x y))
 534
 535     ((integer ratio)
 536      (if (zerop x)
 537          0
 538          (let* ((ny (numerator y))
 539                 (dy (denominator y))
 540                 (gcd (gcd x ny)))
 541            (build-ratio (* (maybe-truncate x gcd) dy)
 542                         (maybe-truncate ny gcd)))))
 543
 544     ((ratio integer)
 545      (let* ((nx (numerator x))
 546             (gcd (gcd nx y)))
 547        (build-ratio (maybe-truncate nx gcd)
 548                     (* (maybe-truncate y gcd) (denominator x)))))))
 549
 550 (defun %negate (n)
 551   (number-dispatch ((n number))
 552     (((foreach fixnum single-float double-float #!+long-float long-float))
 553      (%negate n))
 554     ((bignum)
 555      (negate-bignum n))
 556     ((ratio)
 557      (%make-ratio (- (numerator n)) (denominator n)))
 558     ((complex)
 559      (complex (- (realpart n)) (- (imagpart n))))))
 560 \f
 561 ;;;; TRUNCATE and friends
 562
 563 (defun truncate (number &optional (divisor 1))
 564   #!+sb-doc
 565   "Return number (or number/divisor) as an integer, rounded toward 0.
 566   The second returned value is the remainder."
 567   (macrolet ((truncate-float (rtype)
 568                `(let* ((float-div (coerce divisor ',rtype))
 569                        (res (%unary-truncate (/ number float-div))))
 570                   (values res
 571                           (- number
 572                              (* (coerce res ',rtype) float-div))))))
 573     (number-dispatch ((number real) (divisor real))
 574       ((fixnum fixnum) (truncate number divisor))
 575       (((foreach fixnum bignum) ratio)
 576        (let ((q (truncate (* number (denominator divisor))
 577                           (numerator divisor))))
 578          (values q (- number (* q divisor)))))
 579       ((fixnum bignum)
 580        (values 0 number))
 581       ((ratio (or float rational))
 582        (let ((q (truncate (numerator number)
 583                           (* (denominator number) divisor))))
 584          (values q (- number (* q divisor)))))
 585       ((bignum fixnum)
 586        (bignum-truncate number (make-small-bignum divisor)))
 587       ((bignum bignum)
 588        (bignum-truncate number divisor))
 589
 590       (((foreach single-float double-float #!+long-float long-float)
 591         (or rational single-float))
 592        (if (eql divisor 1)
 593            (let ((res (%unary-truncate number)))
 594              (values res (- number (coerce res '(dispatch-type number)))))
 595            (truncate-float (dispatch-type number))))
 596       #!+long-float
 597       ((long-float (or single-float double-float long-float))
 598        (truncate-float long-float))
 599       #!+long-float
 600       (((foreach double-float single-float) long-float)
 601        (truncate-float long-float))
 602       ((double-float (or single-float double-float))
 603        (truncate-float double-float))
 604       ((single-float double-float)
 605        (truncate-float double-float))
 606       (((foreach fixnum bignum ratio)
 607         (foreach single-float double-float #!+long-float long-float))
 608        (truncate-float (dispatch-type divisor))))))
 609
 610 ;;; Declare these guys inline to let them get optimized a little.
 611 ;;; ROUND and FROUND are not declared inline since they seem too
 612 ;;; obscure and too big to inline-expand by default. Also, this gives
 613 ;;; the compiler a chance to pick off the unary float case. Similarly,
 614 ;;; CEILING and FLOOR are only maybe-inline for now, so that the
 615 ;;; power-of-2 CEILING and FLOOR transforms get a chance.
 616 #!-sb-fluid (declaim (inline rem mod fceiling ffloor ftruncate))
 617 (declaim (maybe-inline ceiling floor))
 618
 619 (defun floor (number &optional (divisor 1))
 620   #!+sb-doc
 621   "Return the greatest integer not greater than number, or number/divisor.
 622   The second returned value is (mod number divisor)."
 623   ;; If the numbers do not divide exactly and the result of
 624   ;; (/ NUMBER DIVISOR) would be negative then decrement the quotient
 625   ;; and augment the remainder by the divisor.
 626   (multiple-value-bind (tru rem) (truncate number divisor)
 627     (if (and (not (zerop rem))
 628              (if (minusp divisor)
 629                  (plusp number)
 630                  (minusp number)))
 631         (values (1- tru) (+ rem divisor))
 632         (values tru rem))))
 633
 634 (defun ceiling (number &optional (divisor 1))
 635   #!+sb-doc
 636   "Return the smallest integer not less than number, or number/divisor.
 637   The second returned value is the remainder."
 638   ;; If the numbers do not divide exactly and the result of
 639   ;; (/ NUMBER DIVISOR) would be positive then increment the quotient
 640   ;; and decrement the remainder by the divisor.
 641   (multiple-value-bind (tru rem) (truncate number divisor)
 642     (if (and (not (zerop rem))
 643              (if (minusp divisor)
 644                  (minusp number)
 645                  (plusp number)))
 646         (values (+ tru 1) (- rem divisor))
 647         (values tru rem))))
 648
 649 (defun round (number &optional (divisor 1))
 650   #!+sb-doc
 651   "Rounds number (or number/divisor) to nearest integer.
 652   The second returned value is the remainder."
 653   (if (eql divisor 1)
 654       (round number)
 655       (multiple-value-bind (tru rem) (truncate number divisor)
 656         (let ((thresh (/ (abs divisor) 2)))
 657           (cond ((or (> rem thresh)
 658                      (and (= rem thresh) (oddp tru)))
 659                  (if (minusp divisor)
 660                      (values (- tru 1) (+ rem divisor))
 661                      (values (+ tru 1) (- rem divisor))))
 662                 ((let ((-thresh (- thresh)))
 663                    (or (< rem -thresh)
 664                        (and (= rem -thresh) (oddp tru))))
 665                  (if (minusp divisor)
 666                      (values (+ tru 1) (- rem divisor))
 667                      (values (- tru 1) (+ rem divisor))))
 668                 (t (values tru rem)))))))
 669
 670 (defun rem (number divisor)
 671   #!+sb-doc
 672   "Return second result of TRUNCATE."
 673   (multiple-value-bind (tru rem) (truncate number divisor)
 674     (declare (ignore tru))
 675     rem))
 676
 677 (defun mod (number divisor)
 678   #!+sb-doc
 679   "Return second result of FLOOR."
 680   (let ((rem (rem number divisor)))
 681     (if (and (not (zerop rem))
 682              (if (minusp divisor)
 683                  (plusp number)
 684                  (minusp number)))
 685         (+ rem divisor)
 686         rem)))
 687
 688 (macrolet ((def (name op doc)
 689              `(defun ,name (number &optional (divisor 1))
 690                 ,doc
 691                 (multiple-value-bind (res rem) (,op number divisor)
 692                   (values (float res (if (floatp rem) rem 1.0)) rem)))))
 693   (def ffloor floor
 694     "Same as FLOOR, but returns first value as a float.")
 695   (def fceiling ceiling
 696     "Same as CEILING, but returns first value as a float." )
 697   (def ftruncate truncate
 698     "Same as TRUNCATE, but returns first value as a float.")
 699   (def fround round
 700     "Same as ROUND, but returns first value as a float."))
 701 \f
 702 ;;;; comparisons
 703
 704 (defun = (number &rest more-numbers)
 705   #!+sb-doc
 706   "Return T if all of its arguments are numerically equal, NIL otherwise."
 707   (do ((nlist more-numbers (cdr nlist)))
 708       ((atom nlist) T)
 709      (declare (list nlist))
 710      (if (not (= (car nlist) number)) (return nil))))
 711
 712 (defun /= (number &rest more-numbers)
 713   #!+sb-doc
 714   "Return T if no two of its arguments are numerically equal, NIL otherwise."
 715   (do* ((head number (car nlist))
 716         (nlist more-numbers (cdr nlist)))
 717        ((atom nlist) t)
 718      (declare (list nlist))
 719      (unless (do* ((nl nlist (cdr nl)))
 720                   ((atom nl) T)
 721                (declare (list nl))
 722                (if (= head (car nl)) (return nil)))
 723        (return nil))))
 724
 725 (defun < (number &rest more-numbers)
 726   #!+sb-doc
 727   "Return T if its arguments are in strictly increasing order, NIL otherwise."
 728   (do* ((n number (car nlist))
 729         (nlist more-numbers (cdr nlist)))
 730        ((atom nlist) t)
 731      (declare (list nlist))
 732      (if (not (< n (car nlist))) (return nil))))
 733
 734 (defun > (number &rest more-numbers)
 735   #!+sb-doc
 736   "Return T if its arguments are in strictly decreasing order, NIL otherwise."
 737   (do* ((n number (car nlist))
 738         (nlist more-numbers (cdr nlist)))
 739        ((atom nlist) t)
 740      (declare (list nlist))
 741      (if (not (> n (car nlist))) (return nil))))
 742
 743 (defun <= (number &rest more-numbers)
 744   #!+sb-doc
 745   "Return T if arguments are in strictly non-decreasing order, NIL otherwise."
 746   (do* ((n number (car nlist))
 747         (nlist more-numbers (cdr nlist)))
 748        ((atom nlist) t)
 749      (declare (list nlist))
 750      (if (not (<= n (car nlist))) (return nil))))
 751
 752 (defun >= (number &rest more-numbers)
 753   #!+sb-doc
 754   "Return T if arguments are in strictly non-increasing order, NIL otherwise."
 755   (do* ((n number (car nlist))
 756         (nlist more-numbers (cdr nlist)))
 757        ((atom nlist) t)
 758      (declare (list nlist))
 759      (if (not (>= n (car nlist))) (return nil))))
 760
 761 (defun max (number &rest more-numbers)
 762   #!+sb-doc
 763   "Return the greatest of its arguments."
 764   (do ((nlist more-numbers (cdr nlist))
 765        (result number))
 766       ((null nlist) (return result))
 767      (declare (list nlist))
 768      (declare (type real number result))
 769      (if (> (car nlist) result) (setq result (car nlist)))))
 770
 771 (defun min (number &rest more-numbers)
 772   #!+sb-doc
 773   "Return the least of its arguments."
 774   (do ((nlist more-numbers (cdr nlist))
 775        (result number))
 776       ((null nlist) (return result))
 777      (declare (list nlist))
 778      (declare (type real number result))
 779      (if (< (car nlist) result) (setq result (car nlist)))))
 780
 781 (eval-when (:compile-toplevel :execute)
 782
 783 ;;; The INFINITE-X-FINITE-Y and INFINITE-Y-FINITE-X args tell us how
 784 ;;; to handle the case when X or Y is a floating-point infinity and
 785 ;;; the other arg is a rational. (Section 12.1.4.1 of the ANSI spec
 786 ;;; says that comparisons are done by converting the float to a
 787 ;;; rational when comparing with a rational, but infinities can't be
 788 ;;; converted to a rational, so we show some initiative and do it this
 789 ;;; way instead.)
 790 (defun basic-compare (op &key infinite-x-finite-y infinite-y-finite-x)
 791   `(((fixnum fixnum) (,op x y))
 792
 793     ((single-float single-float) (,op x y))
 794     #!+long-float
 795     (((foreach single-float double-float long-float) long-float)
 796      (,op (coerce x 'long-float) y))
 797     #!+long-float
 798     ((long-float (foreach single-float double-float))
 799      (,op x (coerce y 'long-float)))
 800     (((foreach single-float double-float) double-float)
 801      (,op (coerce x 'double-float) y))
 802     ((double-float single-float)
 803      (,op x (coerce y 'double-float)))
 804     (((foreach single-float double-float #!+long-float long-float) rational)
 805      (if (eql y 0)
 806          (,op x (coerce 0 '(dispatch-type x)))
 807          (if (float-infinity-p x)
 808              ,infinite-x-finite-y
 809              (,op (rational x) y))))
 810     (((foreach bignum fixnum ratio) float)
 811      (if (float-infinity-p y)
 812          ,infinite-y-finite-x
 813          (,op x (rational y))))))
 814 ) ; EVAL-WHEN
 815
 816 (macrolet ((def-two-arg-</> (name op ratio-arg1 ratio-arg2 &rest cases)
 817              `(defun ,name (x y)
 818                 (number-dispatch ((x real) (y real))
 819                                  (basic-compare
 820                                   ,op
 821                                   :infinite-x-finite-y
 822                                   (,op x (coerce 0 '(dispatch-type x)))
 823                                   :infinite-y-finite-x
 824                                   (,op (coerce 0 '(dispatch-type y)) y))
 825                                  (((foreach fixnum bignum) ratio)
 826                                   (,op x (,ratio-arg2 (numerator y)
 827                                                       (denominator y))))
 828                                  ((ratio integer)
 829                                   (,op (,ratio-arg1 (numerator x)
 830                                                     (denominator x))
 831                                        y))
 832                                  ((ratio ratio)
 833                                   (,op (* (numerator   (truly-the ratio x))
 834                                           (denominator (truly-the ratio y)))
 835                                        (* (numerator   (truly-the ratio y))
 836                                           (denominator (truly-the ratio x)))))
 837                                  ,@cases))))
 838   (def-two-arg-</> two-arg-< < floor ceiling
 839     ((fixnum bignum)
 840      (bignum-plus-p y))
 841     ((bignum fixnum)
 842      (not (bignum-plus-p x)))
 843     ((bignum bignum)
 844      (minusp (bignum-compare x y))))
 845   (def-two-arg-</> two-arg-> > ceiling floor
 846     ((fixnum bignum)
 847      (not (bignum-plus-p y)))
 848     ((bignum fixnum)
 849      (bignum-plus-p x))
 850     ((bignum bignum)
 851      (plusp (bignum-compare x y)))))
 852
 853 (defun two-arg-= (x y)
 854   (number-dispatch ((x number) (y number))
 855     (basic-compare =
 856                    ;; An infinite value is never equal to a finite value.
 857                    :infinite-x-finite-y nil
 858                    :infinite-y-finite-x nil)
 859     ((fixnum (or bignum ratio)) nil)
 860
 861     ((bignum (or fixnum ratio)) nil)
 862     ((bignum bignum)
 863      (zerop (bignum-compare x y)))
 864
 865     ((ratio integer) nil)
 866     ((ratio ratio)
 867      (and (eql (numerator x) (numerator y))
 868           (eql (denominator x) (denominator y))))
 869
 870     ((complex complex)
 871      (and (= (realpart x) (realpart y))
 872           (= (imagpart x) (imagpart y))))
 873     (((foreach fixnum bignum ratio single-float double-float
 874                #!+long-float long-float) complex)
 875      (and (= x (realpart y))
 876           (zerop (imagpart y))))
 877     ((complex (or float rational))
 878      (and (= (realpart x) y)
 879           (zerop (imagpart x))))))
 880
 881 (defun eql (obj1 obj2)
 882   #!+sb-doc
 883   "Return T if OBJ1 and OBJ2 represent the same object, otherwise NIL."
 884   (or (eq obj1 obj2)
 885       (if (or (typep obj2 'fixnum)
 886               (not (typep obj2 'number)))
 887           nil
 888           (macrolet ((foo (&rest stuff)
 889                        `(typecase obj2
 890                           ,@(mapcar (lambda (foo)
 891                                       (let ((type (car foo))
 892                                             (fn (cadr foo)))
 893                                         `(,type
 894                                           (and (typep obj1 ',type)
 895                                                (,fn obj1 obj2)))))
 896                                     stuff))))
 897             (foo
 898               (single-float eql)
 899               (double-float eql)
 900               #!+long-float
 901               (long-float eql)
 902               (bignum
 903                (lambda (x y)
 904                  (zerop (bignum-compare x y))))
 905               (ratio
 906                (lambda (x y)
 907                  (and (eql (numerator x) (numerator y))
 908                       (eql (denominator x) (denominator y)))))
 909               (complex
 910                (lambda (x y)
 911                  (and (eql (realpart x) (realpart y))
 912                       (eql (imagpart x) (imagpart y))))))))))
 913 \f
 914 ;;;; logicals
 915
 916 (defun logior (&rest integers)
 917   #!+sb-doc
 918   "Return the bit-wise or of its arguments. Args must be integers."
 919   (declare (list integers))
 920   (if integers
 921       (do ((result (pop integers) (logior result (pop integers))))
 922           ((null integers) result)
 923         (declare (integer result)))
 924       0))
 925
 926 (defun logxor (&rest integers)
 927   #!+sb-doc
 928   "Return the bit-wise exclusive or of its arguments. Args must be integers."
 929   (declare (list integers))
 930   (if integers
 931       (do ((result (pop integers) (logxor result (pop integers))))
 932           ((null integers) result)
 933         (declare (integer result)))
 934       0))
 935
 936 (defun logand (&rest integers)
 937   #!+sb-doc
 938   "Return the bit-wise and of its arguments. Args must be integers."
 939   (declare (list integers))
 940   (if integers
 941       (do ((result (pop integers) (logand result (pop integers))))
 942           ((null integers) result)
 943         (declare (integer result)))
 944       -1))
 945
 946 (defun logeqv (&rest integers)
 947   #!+sb-doc
 948   "Return the bit-wise equivalence of its arguments. Args must be integers."
 949   (declare (list integers))
 950   (if integers
 951       (do ((result (pop integers) (logeqv result (pop integers))))
 952           ((null integers) result)
 953         (declare (integer result)))
 954       -1))
 955
 956 (defun lognand (integer1 integer2)
 957   #!+sb-doc
 958   "Return the complement of the logical AND of integer1 and integer2."
 959   (lognand integer1 integer2))
 960
 961 (defun lognor (integer1 integer2)
 962   #!+sb-doc
 963   "Return the complement of the logical OR of integer1 and integer2."
 964   (lognor integer1 integer2))
 965
 966 (defun logandc1 (integer1 integer2)
 967   #!+sb-doc
 968   "Return the logical AND of (LOGNOT integer1) and integer2."
 969   (logandc1 integer1 integer2))
 970
 971 (defun logandc2 (integer1 integer2)
 972   #!+sb-doc
 973   "Return the logical AND of integer1 and (LOGNOT integer2)."
 974   (logandc2 integer1 integer2))
 975
 976 (defun logorc1 (integer1 integer2)
 977   #!+sb-doc
 978   "Return the logical OR of (LOGNOT integer1) and integer2."
 979   (logorc1 integer1 integer2))
 980
 981 (defun logorc2 (integer1 integer2)
 982   #!+sb-doc
 983   "Return the logical OR of integer1 and (LOGNOT integer2)."
 984   (logorc2 integer1 integer2))
 985
 986 (defun lognot (number)
 987   #!+sb-doc
 988   "Return the bit-wise logical not of integer."
 989   (etypecase number
 990     (fixnum (lognot (truly-the fixnum number)))
 991     (bignum (bignum-logical-not number))))
 992
 993 (macrolet ((def (name op big-op)
 994              `(defun ,name (x y)
 995                (number-dispatch ((x integer) (y integer))
 996                  (bignum-cross-fixnum ,op ,big-op)))))
 997   (def two-arg-and logand bignum-logical-and)
 998   (def two-arg-ior logior bignum-logical-ior)
 999   (def two-arg-xor logxor bignum-logical-xor))
1000
1001 (defun logcount (integer)
1002   #!+sb-doc
1003   "Count the number of 1 bits if INTEGER is positive, and the number of 0 bits
1004   if INTEGER is negative."
1005   (etypecase integer
1006     (fixnum
1007      (logcount (truly-the (integer 0 #.(max most-positive-fixnum
1008                                             (lognot most-negative-fixnum)))
1009                           (if (minusp (truly-the fixnum integer))
1010                               (lognot (truly-the fixnum integer))
1011                               integer))))
1012     (bignum
1013      (bignum-logcount integer))))
1014
1015 (defun logtest (integer1 integer2)
1016   #!+sb-doc
1017   "Predicate which returns T if logand of integer1 and integer2 is not zero."
1018   (logtest integer1 integer2))
1019
1020 (defun logbitp (index integer)
1021   #!+sb-doc
1022   "Predicate returns T if bit index of integer is a 1."
1023   (logbitp index integer))
1024
1025 (defun ash (integer count)
1026   #!+sb-doc
1027   "Shifts integer left by count places preserving sign. - count shifts right."
1028   (declare (integer integer count))
1029   (etypecase integer
1030     (fixnum
1031      (cond ((zerop integer)
1032             0)
1033            ((fixnump count)
1034             (let ((length (integer-length (truly-the fixnum integer)))
1035                   (count (truly-the fixnum count)))
1036               (declare (fixnum length count))
1037               (cond ((and (plusp count)
1038                           (> (+ length count)
1039                              (integer-length most-positive-fixnum)))
1040                      (bignum-ashift-left (make-small-bignum integer) count))
1041                     (t
1042                      (truly-the fixnum
1043                                 (ash (truly-the fixnum integer) count))))))
1044            ((minusp count)
1045             (if (minusp integer) -1 0))
1046            (t
1047             (bignum-ashift-left (make-small-bignum integer) count))))
1048     (bignum
1049      (if (plusp count)
1050          (bignum-ashift-left integer count)
1051          (bignum-ashift-right integer (- count))))))
1052
1053 (defun integer-length (integer)
1054   #!+sb-doc
1055   "Return the number of significant bits in the absolute value of integer."
1056   (etypecase integer
1057     (fixnum
1058      (integer-length (truly-the fixnum integer)))
1059     (bignum
1060      (bignum-integer-length integer))))
1061 \f
1062 ;;;; BYTE, bytespecs, and related operations
1063
1064 (defun byte (size position)
1065   #!+sb-doc
1066   "Return a byte specifier which may be used by other byte functions
1067   (e.g. LDB)."
1068   (byte size position))
1069
1070 (defun byte-size (bytespec)
1071   #!+sb-doc
1072   "Return the size part of the byte specifier bytespec."
1073   (byte-size bytespec))
1074
1075 (defun byte-position (bytespec)
1076   #!+sb-doc
1077   "Return the position part of the byte specifier bytespec."
1078   (byte-position bytespec))
1079
1080 (defun ldb (bytespec integer)
1081   #!+sb-doc
1082   "Extract the specified byte from integer, and right justify result."
1083   (ldb bytespec integer))
1084
1085 (defun ldb-test (bytespec integer)
1086   #!+sb-doc
1087   "Return T if any of the specified bits in integer are 1's."
1088   (ldb-test bytespec integer))
1089
1090 (defun mask-field (bytespec integer)
1091   #!+sb-doc
1092   "Extract the specified byte from integer,  but do not right justify result."
1093   (mask-field bytespec integer))
1094
1095 (defun dpb (newbyte bytespec integer)
1096   #!+sb-doc
1097   "Return new integer with newbyte in specified position, newbyte is right justified."
1098   (dpb newbyte bytespec integer))
1099
1100 (defun deposit-field (newbyte bytespec integer)
1101   #!+sb-doc
1102   "Return new integer with newbyte in specified position, newbyte is not right justified."
1103   (deposit-field newbyte bytespec integer))
1104
1105 (defun %ldb (size posn integer)
1106   (logand (ash integer (- posn))
1107           (1- (ash 1 size))))
1108
1109 (defun %mask-field (size posn integer)
1110   (logand integer (ash (1- (ash 1 size)) posn)))
1111
1112 (defun %dpb (newbyte size posn integer)
1113   (let ((mask (1- (ash 1 size))))
1114     (logior (logand integer (lognot (ash mask posn)))
1115             (ash (logand newbyte mask) posn))))
1116
1117 (defun %deposit-field (newbyte size posn integer)
1118   (let ((mask (ash (ldb (byte size 0) -1) posn)))
1119     (logior (logand newbyte mask)
1120             (logand integer (lognot mask)))))
1121 \f
1122 ;;;; BOOLE
1123
1124 ;;; The boole function dispaches to any logic operation depending on
1125 ;;;     the value of a variable. Presently, legal selector values are [0..15].
1126 ;;;     boole is open coded for calls with a constant selector. or with calls
1127 ;;;     using any of the constants declared below.
1128
1129 (defconstant boole-clr 0
1130   #!+sb-doc
1131   "Boole function op, makes BOOLE return 0.")
1132
1133 (defconstant boole-set 1
1134   #!+sb-doc
1135   "Boole function op, makes BOOLE return -1.")
1136
1137 (defconstant boole-1   2
1138   #!+sb-doc
1139   "Boole function op, makes BOOLE return integer1.")
1140
1141 (defconstant boole-2   3
1142   #!+sb-doc
1143   "Boole function op, makes BOOLE return integer2.")
1144
1145 (defconstant boole-c1  4
1146   #!+sb-doc
1147   "Boole function op, makes BOOLE return complement of integer1.")
1148
1149 (defconstant boole-c2  5
1150   #!+sb-doc
1151   "Boole function op, makes BOOLE return complement of integer2.")
1152
1153 (defconstant boole-and 6
1154   #!+sb-doc
1155   "Boole function op, makes BOOLE return logand of integer1 and integer2.")
1156
1157 (defconstant boole-ior 7
1158   #!+sb-doc
1159   "Boole function op, makes BOOLE return logior of integer1 and integer2.")
1160
1161 (defconstant boole-xor 8
1162   #!+sb-doc
1163   "Boole function op, makes BOOLE return logxor of integer1 and integer2.")
1164
1165 (defconstant boole-eqv 9
1166   #!+sb-doc
1167   "Boole function op, makes BOOLE return logeqv of integer1 and integer2.")
1168
1169 (defconstant boole-nand  10
1170   #!+sb-doc
1171   "Boole function op, makes BOOLE return log nand of integer1 and integer2.")
1172
1173 (defconstant boole-nor   11
1174   #!+sb-doc
1175   "Boole function op, makes BOOLE return lognor of integer1 and integer2.")
1176
1177 (defconstant boole-andc1 12
1178   #!+sb-doc
1179   "Boole function op, makes BOOLE return logandc1 of integer1 and integer2.")
1180
1181 (defconstant boole-andc2 13
1182   #!+sb-doc
1183   "Boole function op, makes BOOLE return logandc2 of integer1 and integer2.")
1184
1185 (defconstant boole-orc1  14
1186   #!+sb-doc
1187   "Boole function op, makes BOOLE return logorc1 of integer1 and integer2.")
1188
1189 (defconstant boole-orc2  15
1190   #!+sb-doc
1191   "Boole function op, makes BOOLE return logorc2 of integer1 and integer2.")
1192
1193 (defun boole (op integer1 integer2)
1194   #!+sb-doc
1195   "Bit-wise boolean function on two integers. Function chosen by OP:
1196         0       BOOLE-CLR
1197         1       BOOLE-SET
1198         2       BOOLE-1
1199         3       BOOLE-2
1200         4       BOOLE-C1
1201         5       BOOLE-C2
1202         6       BOOLE-AND
1203         7       BOOLE-IOR
1204         8       BOOLE-XOR
1205         9       BOOLE-EQV
1206         10      BOOLE-NAND
1207         11      BOOLE-NOR
1208         12      BOOLE-ANDC1
1209         13      BOOLE-ANDC2
1210         14      BOOLE-ORC1
1211         15      BOOLE-ORC2"
1212   (case op
1213     (0 (boole 0 integer1 integer2))
1214     (1 (boole 1 integer1 integer2))
1215     (2 (boole 2 integer1 integer2))
1216     (3 (boole 3 integer1 integer2))
1217     (4 (boole 4 integer1 integer2))
1218     (5 (boole 5 integer1 integer2))
1219     (6 (boole 6 integer1 integer2))
1220     (7 (boole 7 integer1 integer2))
1221     (8 (boole 8 integer1 integer2))
1222     (9 (boole 9 integer1 integer2))
1223     (10 (boole 10 integer1 integer2))
1224     (11 (boole 11 integer1 integer2))
1225     (12 (boole 12 integer1 integer2))
1226     (13 (boole 13 integer1 integer2))
1227     (14 (boole 14 integer1 integer2))
1228     (15 (boole 15 integer1 integer2))
1229     (t (error 'type-error :datum op :expected-type '(mod 16)))))
1230 \f
1231 ;;;; GCD and LCM
1232
1233 (defun gcd (&rest numbers)
1234   #!+sb-doc
1235   "Return the greatest common divisor of the arguments, which must be
1236   integers. Gcd with no arguments is defined to be 0."
1237   (cond ((null numbers) 0)
1238         ((null (cdr numbers)) (abs (the integer (car numbers))))
1239         (t
1240          (do ((gcd (the integer (car numbers))
1241                    (gcd gcd (the integer (car rest))))
1242               (rest (cdr numbers) (cdr rest)))
1243              ((null rest) gcd)
1244            (declare (integer gcd)
1245                     (list rest))))))
1246
1247 (defun lcm (&rest numbers)
1248   #!+sb-doc
1249   "Return the least common multiple of one or more integers. LCM of no
1250   arguments is defined to be 1."
1251   (cond ((null numbers) 1)
1252         ((null (cdr numbers)) (abs (the integer (car numbers))))
1253         (t
1254          (do ((lcm (the integer (car numbers))
1255                    (lcm lcm (the integer (car rest))))
1256               (rest (cdr numbers) (cdr rest)))
1257              ((null rest) lcm)
1258            (declare (integer lcm) (list rest))))))
1259
1260 (defun two-arg-lcm (n m)
1261   (declare (integer n m))
1262   (* (truncate (max n m) (gcd n m)) (min n m)))
1263
1264 ;;; Do the GCD of two integer arguments. With fixnum arguments, we use the
1265 ;;; binary GCD algorithm from Knuth's seminumerical algorithms (slightly
1266 ;;; structurified), otherwise we call BIGNUM-GCD. We pick off the special case
1267 ;;; of 0 before the dispatch so that the bignum code doesn't have to worry
1268 ;;; about "small bignum" zeros.
1269 (defun two-arg-gcd (u v)
1270   (cond ((eql u 0) v)
1271         ((eql v 0) u)
1272         (t
1273          (number-dispatch ((u integer) (v integer))
1274            ((fixnum fixnum)
1275             (locally
1276               (declare (optimize (speed 3) (safety 0)))
1277               (do ((k 0 (1+ k))
1278                    (u (abs u) (ash u -1))
1279                    (v (abs v) (ash v -1)))
1280                   ((oddp (logior u v))
1281                    (do ((temp (if (oddp u) (- v) (ash u -1))
1282                               (ash temp -1)))
1283                        (nil)
1284                      (declare (fixnum temp))
1285                      (when (oddp temp)
1286                        (if (plusp temp)
1287                            (setq u temp)
1288                            (setq v (- temp)))
1289                        (setq temp (- u v))
1290                        (when (zerop temp)
1291                          (let ((res (ash u k)))
1292                            (declare (type (signed-byte 31) res)
1293                                     (optimize (inhibit-warnings 3)))
1294                            (return res))))))
1295                 (declare (type (mod 30) k)
1296                          (type (signed-byte 31) u v)))))
1297            ((bignum bignum)
1298             (bignum-gcd u v))
1299            ((bignum fixnum)
1300             (bignum-gcd u (make-small-bignum v)))
1301            ((fixnum bignum)
1302             (bignum-gcd (make-small-bignum u) v))))))
1303 \f
1304 ;;; From discussion on comp.lang.lisp and Akira Kurihara.
1305 (defun isqrt (n)
1306   #!+sb-doc
1307   "Return the root of the nearest integer less than n which is a perfect
1308    square."
1309   (declare (type unsigned-byte n) (values unsigned-byte))
1310   ;; Theoretically (> n 7), i.e., n-len-quarter > 0.
1311   (if (and (fixnump n) (<= n 24))
1312       (cond ((> n 15) 4)
1313             ((> n  8) 3)
1314             ((> n  3) 2)
1315             ((> n  0) 1)
1316             (t 0))
1317       (let* ((n-len-quarter (ash (integer-length n) -2))
1318              (n-half (ash n (- (ash n-len-quarter 1))))
1319              (n-half-isqrt (isqrt n-half))
1320              (init-value (ash (1+ n-half-isqrt) n-len-quarter)))
1321         (loop
1322           (let ((iterated-value
1323                  (ash (+ init-value (truncate n init-value)) -1)))
1324             (unless (< iterated-value init-value)
1325               (return init-value))
1326             (setq init-value iterated-value))))))
1327 \f
1328 ;;;; miscellaneous number predicates
1329
1330 (macrolet ((def (name doc)
1331              `(defun ,name (number) ,doc (,name number))))
1332   (def zerop "Is this number zero?")
1333   (def plusp "Is this real number strictly positive?")
1334   (def minusp "Is this real number strictly negative?")
1335   (def oddp "Is this integer odd?")
1336   (def evenp "Is this integer even?"))