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