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