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