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