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