1 ;;;; This file contains the definitions of most number functions.
3 ;;;; This software is part of the SBCL system. See the README file for
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.
12 (in-package "SB!KERNEL")
14 ;;;; the NUMBER-DISPATCH macro
16 (eval-when (:compile-toplevel
:load-toplevel
:execute
)
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
)
24 (unless (null types
) (error "More types than vars."))
26 (error "Duplicate case: ~S." body
))
28 (sublis var-types body
:test
#'equal
)))
30 (error "More vars than types."))
32 (flet ((frob (var type
)
33 (parse-number-dispatch
35 (or (assoc type
(cdr result
) :test
#'equal
)
36 (car (setf (cdr result
)
37 (acons type nil
(cdr result
)))))
39 (acons `(dispatch-type ,var
) type var-types
)
41 (let ((type (first types
))
43 (if (and (consp type
) (eq (first type
) 'foreach
))
44 (dolist (type (rest type
))
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
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
*)))
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
)
67 (let ((var (first vars
))
68 (cases (sort cases
#'type-test-order
:key
#'car
)))
70 ,@(mapcar (lambda (case)
72 ,@(generate-number-dispatch (rest vars
)
76 (t (go ,(first error-tags
))))))
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.
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
))
100 (if (symbolp (first case
))
101 (let ((cases (apply (symbol-function (first case
)) (rest case
))))
103 (parse-number-dispatch vars res
(first case
) nil
(rest case
))))
104 (parse-number-dispatch vars res
(first case
) nil
(rest case
))))
108 (dolist (spec var-specs
)
109 (let ((var (first spec
))
114 (errors `(return-from
116 (error 'simple-type-error
:datum
,var
117 :expected-type
',type
119 "~@<Argument ~A is not a ~S: ~2I~_~S~:>"
121 (list ',var
',type
,var
))))))
126 ,@(generate-number-dispatch vars
(error-tags)
130 ;;;; binary operation dispatching utilities
132 (eval-when (:compile-toplevel
:execute
)
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
))))
144 (((foreach single-float double-float long-float
) long-float
)
145 (,op
(coerce ,x
'long-float
) ,y
))
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
)))))
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
))
158 (,big-op
(make-small-bignum x
) y
))
160 (,big-op x
(make-small-bignum y
)))
166 ;;;; canonicalization utilities
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
)
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
)))
186 (%make-complex realpart imagpart
)))))
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
)
195 (values (- num
) (- den
))
199 (error 'division-by-zero
200 :operands
(list num den
)
201 :operation
'build-ratio
))
203 (t (%make-ratio num den
)))))
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
)
214 (defun upgraded-complex-part-type (spec &optional environment
)
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
)
223 ((subtypep spec
'double-float
)
226 ((subtypep spec
'long-float
)
228 ((subtypep spec
'rational
)
233 (defun complex (realpart &optional
(imagpart 0))
235 "Return a complex number with the specified real and imaginary components."
236 (flet ((%%make-complex
(realpart imagpart
)
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
)))
251 (%make-complex realpart imagpart
)))))
252 (number-dispatch ((realpart real
) (imagpart real
))
254 (canonical-complex realpart imagpart
))
255 (float-contagion %%make-complex realpart imagpart
(rational)))))
257 (defun realpart (number)
259 "Extract the real part of a number."
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
)))
269 (sb!kernel
:%realpart number
))
273 (defun imagpart (number)
275 "Extract the imaginary part of a number."
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
)))
285 (sb!kernel
:%imagpart number
))
291 (defun conjugate (number)
293 "Return the complex conjugate of NUMBER. For non-complex numbers, this is
295 (if (complexp number
)
296 (complex (realpart number
) (- (imagpart number
)))
299 (defun signum (number)
301 "If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))."
304 (if (rationalp number
)
305 (if (plusp number
) 1 -
1)
306 (/ number
(abs number
)))))
310 (defun numerator (number)
312 "Return the numerator of NUMBER, which must be rational."
315 (defun denominator (number)
317 "Return the denominator of NUMBER, which must be rational."
318 (denominator number
))
320 ;;;; arithmetic operations
322 (macrolet ((define-arith (op init doc
)
323 #!-sb-doc
(declare (ignore doc
))
324 `(defun ,op
(&rest args
)
326 (if (null args
) ,init
327 (do ((args (cdr args
) (cdr args
))
328 (result (car args
) (,op result
(car args
))))
330 ;; to signal TYPE-ERROR when exactly 1 arg of wrong type:
331 (declare (type number result
)))))))
333 "Return the sum of its arguments. With no args, returns 0.")
335 "Return the product of its arguments. With no args, returns 1."))
337 (defun - (number &rest more-numbers
)
339 "Subtract the second and all subsequent arguments from the first;
340 or with one argument, negate the first argument."
342 (do ((nlist more-numbers
(cdr nlist
))
344 ((atom nlist
) result
)
345 (declare (list nlist
))
346 (setq result
(- result
(car nlist
))))
349 (defun / (number &rest more-numbers
)
351 "Divide the first argument by each of the following arguments, in turn.
352 With one argument, return reciprocal."
354 (do ((nlist more-numbers
(cdr nlist
))
356 ((atom nlist
) result
)
357 (declare (list nlist
))
358 (setq result
(/ result
(car nlist
))))
371 (eval-when (:compile-toplevel
)
373 (sb!xc
:defmacro two-arg-
+/-
(name op big-op
)
375 (number-dispatch ((x number
) (y number
))
376 (bignum-cross-fixnum ,op
,big-op
)
377 (float-contagion ,op x y
)
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
)))
388 (((foreach fixnum bignum
) ratio
)
389 (let* ((dy (denominator y
))
390 (n (,op
(* x dy
) (numerator y
))))
393 (let* ((dx (denominator x
))
394 (n (,op
(numerator x
) (* y dx
))))
397 (let* ((nx (numerator x
))
403 (%make-ratio
(,op
(* nx dy
) (* dx ny
)) (* dx dy
))
404 (let* ((t1 (,op
(* nx
(truncate dy g1
)) (* (truncate dx g1
) ny
)))
406 (t2 (truncate dx g1
)))
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
))))))))))))
417 (two-arg-+/- two-arg-
+ + add-bignums
)
418 (two-arg-+/- two-arg-- - subtract-bignum
)
420 (defun two-arg-* (x y
)
421 (flet ((integer*ratio
(x y
)
423 (let* ((ny (numerator y
))
427 (%make-ratio
(* x ny
) dy
)
428 (let ((nn (* (truncate x gcd
) ny
))
429 (nd (truncate dy gcd
)))
432 (%make-ratio nn nd
)))))))
434 (canonical-complex (* (realpart x
) y
) (* (imagpart x
) y
))))
435 (number-dispatch ((x number
) (y number
))
436 (float-contagion * x y
)
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
))
444 (let* ((rx (realpart x
))
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
)
453 ((complex (or rational float
))
456 (((foreach bignum fixnum
) ratio
) (integer*ratio x y
))
457 ((ratio integer
) (integer*ratio y x
))
459 (let* ((nx (numerator x
))
465 (build-ratio (* (maybe-truncate nx g1
)
466 (maybe-truncate ny g2
))
467 (* (maybe-truncate dx g2
)
468 (maybe-truncate dy g1
))))))))
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
)
479 (let ((gcd (gcd x y
)))
480 (declare (fixnum gcd
))
483 (build-ratio (truncate x gcd
) (truncate y gcd
))))))
484 (let ((gcd (gcd x y
)))
487 (build-ratio (truncate x gcd
) (truncate y gcd
))))))
489 (defun two-arg-/ (x y
)
490 (number-dispatch ((x number
) (y number
))
491 (float-contagion / x y
(ratio integer
))
494 (let* ((rx (realpart x
))
498 (if (> (abs ry
) (abs iy
))
500 (dn (* ry
(+ 1 (* r r
)))))
501 (canonical-complex (/ (+ rx
(* ix r
)) dn
)
502 (/ (- ix
(* rx r
)) dn
)))
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
))
510 (if (> (abs ry
) (abs iy
))
512 (dn (* ry
(+ 1 (* r r
)))))
513 (canonical-complex (/ x dn
)
516 (dn (* iy
(+ 1 (* r r
)))))
517 (canonical-complex (/ (* x r
) dn
)
519 ((complex (or rational float
))
520 (canonical-complex (/ (realpart x
) y
)
524 (let* ((nx (numerator x
))
530 (build-ratio (* (maybe-truncate nx g1
) (maybe-truncate dy g2
))
531 (* (maybe-truncate dx g2
) (maybe-truncate ny g1
)))))
534 (integer-/-integer x y
))
539 (let* ((ny (numerator y
))
542 (build-ratio (* (maybe-truncate x gcd
) dy
)
543 (maybe-truncate ny gcd
)))))
546 (let* ((nx (numerator x
))
548 (build-ratio (maybe-truncate nx gcd
)
549 (* (maybe-truncate y gcd
) (denominator x
)))))))
552 (number-dispatch ((n number
))
553 (((foreach fixnum single-float double-float
#!+long-float long-float
))
558 (%make-ratio
(- (numerator n
)) (denominator n
)))
560 (complex (- (realpart n
)) (- (imagpart n
))))))
562 ;;;; TRUNCATE and friends
564 (defun truncate (number &optional
(divisor 1))
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
))))
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
)))))
581 (bignum-truncate (make-small-bignum number
) divisor
))
582 ((ratio (or float rational
))
583 (let ((q (truncate (numerator number
)
584 (* (denominator number
) divisor
))))
585 (values q
(- number
(* q divisor
)))))
587 (bignum-truncate number
(make-small-bignum divisor
)))
589 (bignum-truncate number divisor
))
591 (((foreach single-float double-float
#!+long-float long-float
)
592 (or rational single-float
))
594 (let ((res (%unary-truncate number
)))
595 (values res
(- number
(coerce res
'(dispatch-type number
)))))
596 (truncate-float (dispatch-type number
))))
598 ((long-float (or single-float double-float long-float
))
599 (truncate-float 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
))))))
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
))
620 (defun floor (number &optional
(divisor 1))
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
))
632 (values (1- tru
) (+ rem divisor
))
635 (defun ceiling (number &optional
(divisor 1))
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
))
647 (values (+ tru
1) (- rem divisor
))
650 (defun round (number &optional
(divisor 1))
652 "Rounds number (or number/divisor) to nearest integer.
653 The second returned value is the remainder."
656 (multiple-value-bind (tru rem
) (truncate number divisor
)
659 (let ((thresh (/ (abs divisor
) 2)))
660 (cond ((or (> rem thresh
)
661 (and (= rem thresh
) (oddp tru
)))
663 (values (- tru
1) (+ rem divisor
))
664 (values (+ tru
1) (- rem divisor
))))
665 ((let ((-thresh (- thresh
)))
667 (and (= rem -thresh
) (oddp tru
))))
669 (values (+ tru
1) (- rem divisor
))
670 (values (- tru
1) (+ rem divisor
))))
671 (t (values tru rem
))))))))
673 (defun rem (number divisor
)
675 "Return second result of TRUNCATE."
676 (multiple-value-bind (tru rem
) (truncate number divisor
)
677 (declare (ignore tru
))
680 (defun mod (number divisor
)
682 "Return second result of FLOOR."
683 (let ((rem (rem number divisor
)))
684 (if (and (not (zerop rem
))
691 (defmacro !define-float-rounding-function
(name op doc
)
692 `(defun ,name
(number &optional
(divisor 1))
694 (multiple-value-bind (res rem
) (,op number divisor
)
695 (values (float res
(if (floatp rem
) rem
1.0)) rem
))))
697 (!define-float-rounding-function ffloor floor
698 "Same as FLOOR, but returns first value as a float.")
699 (!define-float-rounding-function fceiling ceiling
700 "Same as CEILING, but returns first value as a float." )
701 (!define-float-rounding-function ftruncate truncate
702 "Same as TRUNCATE, but returns first value as a float.")
703 (!define-float-rounding-function fround round
704 "Same as ROUND, but returns first value as a float.")
708 (defun = (number &rest more-numbers
)
710 "Return T if all of its arguments are numerically equal, NIL otherwise."
711 (do ((nlist more-numbers
(cdr nlist
)))
713 (declare (list nlist
))
714 (if (not (= (car nlist
) number
)) (return nil
))))
716 (defun /= (number &rest more-numbers
)
718 "Return T if no two of its arguments are numerically equal, NIL otherwise."
719 (do* ((head number
(car nlist
))
720 (nlist more-numbers
(cdr nlist
)))
722 (declare (list nlist
))
723 (unless (do* ((nl nlist
(cdr nl
)))
726 (if (= head
(car nl
)) (return nil
)))
729 (defun < (number &rest more-numbers
)
731 "Return T if its arguments are in strictly increasing order, NIL otherwise."
732 (do* ((n number
(car nlist
))
733 (nlist more-numbers
(cdr nlist
)))
735 (declare (list nlist
))
736 (if (not (< n
(car nlist
))) (return nil
))))
738 (defun > (number &rest more-numbers
)
740 "Return T if its arguments are in strictly decreasing order, NIL otherwise."
741 (do* ((n number
(car nlist
))
742 (nlist more-numbers
(cdr nlist
)))
744 (declare (list nlist
))
745 (if (not (> n
(car nlist
))) (return nil
))))
747 (defun <= (number &rest more-numbers
)
749 "Return T if arguments are in strictly non-decreasing order, NIL otherwise."
750 (do* ((n number
(car nlist
))
751 (nlist more-numbers
(cdr nlist
)))
753 (declare (list nlist
))
754 (if (not (<= n
(car nlist
))) (return nil
))))
756 (defun >= (number &rest more-numbers
)
758 "Return T if arguments are in strictly non-increasing order, NIL otherwise."
759 (do* ((n number
(car nlist
))
760 (nlist more-numbers
(cdr nlist
)))
762 (declare (list nlist
))
763 (if (not (>= n
(car nlist
))) (return nil
))))
765 (defun max (number &rest more-numbers
)
767 "Return the greatest of its arguments."
768 (do ((nlist more-numbers
(cdr nlist
))
770 ((null nlist
) (return result
))
771 (declare (list nlist
))
772 (declare (type real number result
))
773 (if (> (car nlist
) result
) (setq result
(car nlist
)))))
775 (defun min (number &rest more-numbers
)
777 "Return the least of its arguments."
778 (do ((nlist more-numbers
(cdr nlist
))
780 ((null nlist
) (return result
))
781 (declare (list nlist
))
782 (declare (type real number result
))
783 (if (< (car nlist
) result
) (setq result
(car nlist
)))))
785 (eval-when (:compile-toplevel
:execute
)
787 ;;; The INFINITE-X-FINITE-Y and INFINITE-Y-FINITE-X args tell us how
788 ;;; to handle the case when X or Y is a floating-point infinity and
789 ;;; the other arg is a rational. (Section 12.1.4.1 of the ANSI spec
790 ;;; says that comparisons are done by converting the float to a
791 ;;; rational when comparing with a rational, but infinities can't be
792 ;;; converted to a rational, so we show some initiative and do it this
794 (defun basic-compare (op &key infinite-x-finite-y infinite-y-finite-x
)
795 `(((fixnum fixnum
) (,op x y
))
797 ((single-float single-float
) (,op x y
))
799 (((foreach single-float double-float long-float
) long-float
)
800 (,op
(coerce x
'long-float
) y
))
802 ((long-float (foreach single-float double-float
))
803 (,op x
(coerce y
'long-float
)))
804 (((foreach single-float double-float
) double-float
)
805 (,op
(coerce x
'double-float
) y
))
806 ((double-float single-float
)
807 (,op x
(coerce y
'double-float
)))
808 (((foreach single-float double-float
#!+long-float long-float
) rational
)
810 (,op x
(coerce 0 '(dispatch-type x
)))
811 (if (float-infinity-p x
)
813 (,op
(rational x
) y
))))
814 (((foreach bignum fixnum ratio
) float
)
815 (if (float-infinity-p y
)
817 (,op x
(rational y
))))))
820 (macrolet ((def-two-arg-</> (name op ratio-arg1 ratio-arg2
&rest cases
)
822 (number-dispatch ((x real
) (y real
))
826 (,op x
(coerce 0 '(dispatch-type x
)))
828 (,op
(coerce 0 '(dispatch-type y
)) y
))
829 (((foreach fixnum bignum
) ratio
)
830 (,op x
(,ratio-arg2
(numerator y
)
833 (,op
(,ratio-arg1
(numerator x
)
837 (,op
(* (numerator (truly-the ratio x
))
838 (denominator (truly-the ratio y
)))
839 (* (numerator (truly-the ratio y
))
840 (denominator (truly-the ratio x
)))))
842 (def-two-arg-</> two-arg-
< < floor ceiling
846 (not (bignum-plus-p x
)))
848 (minusp (bignum-compare x y
))))
849 (def-two-arg-</> two-arg-
> > ceiling floor
851 (not (bignum-plus-p y
)))
855 (plusp (bignum-compare x y
)))))
857 (defun two-arg-= (x y
)
858 (number-dispatch ((x number
) (y number
))
860 ;; An infinite value is never equal to a finite value.
861 :infinite-x-finite-y nil
862 :infinite-y-finite-x nil
)
863 ((fixnum (or bignum ratio
)) nil
)
865 ((bignum (or fixnum ratio
)) nil
)
867 (zerop (bignum-compare x y
)))
869 ((ratio integer
) nil
)
871 (and (eql (numerator x
) (numerator y
))
872 (eql (denominator x
) (denominator y
))))
875 (and (= (realpart x
) (realpart y
))
876 (= (imagpart x
) (imagpart y
))))
877 (((foreach fixnum bignum ratio single-float double-float
878 #!+long-float long-float
) complex
)
879 (and (= x
(realpart y
))
880 (zerop (imagpart y
))))
881 ((complex (or float rational
))
882 (and (= (realpart x
) y
)
883 (zerop (imagpart x
))))))
885 (defun eql (obj1 obj2
)
887 "Return T if OBJ1 and OBJ2 represent the same object, otherwise NIL."
889 (if (or (typep obj2
'fixnum
)
890 (not (typep obj2
'number
)))
892 (macrolet ((foo (&rest stuff
)
894 ,@(mapcar (lambda (foo)
895 (let ((type (car foo
))
898 (and (typep obj1
',type
)
908 (zerop (bignum-compare x y
))))
911 (and (eql (numerator x
) (numerator y
))
912 (eql (denominator x
) (denominator y
)))))
915 (and (eql (realpart x
) (realpart y
))
916 (eql (imagpart x
) (imagpart y
))))))))))
920 (defun logior (&rest integers
)
922 "Return the bit-wise or of its arguments. Args must be integers."
923 (declare (list integers
))
925 (do ((result (pop integers
) (logior result
(pop integers
))))
926 ((null integers
) result
)
927 (declare (integer result
)))
930 (defun logxor (&rest integers
)
932 "Return the bit-wise exclusive or of its arguments. Args must be integers."
933 (declare (list integers
))
935 (do ((result (pop integers
) (logxor result
(pop integers
))))
936 ((null integers
) result
)
937 (declare (integer result
)))
940 (defun logand (&rest integers
)
942 "Return the bit-wise and of its arguments. Args must be integers."
943 (declare (list integers
))
945 (do ((result (pop integers
) (logand result
(pop integers
))))
946 ((null integers
) result
)
947 (declare (integer result
)))
950 (defun logeqv (&rest integers
)
952 "Return the bit-wise equivalence of its arguments. Args must be integers."
953 (declare (list integers
))
955 (do ((result (pop integers
) (logeqv result
(pop integers
))))
956 ((null integers
) result
)
957 (declare (integer result
)))
960 (defun lognot (number)
962 "Return the bit-wise logical not of integer."
964 (fixnum (lognot (truly-the fixnum number
)))
965 (bignum (bignum-logical-not number
))))
967 (macrolet ((def (name op big-op
&optional doc
)
968 `(defun ,name
(integer1 integer2
)
973 (number-dispatch ((x integer
) (y integer
))
974 (bignum-cross-fixnum ,op
,big-op
))))))
975 (def two-arg-and logand bignum-logical-and
)
976 (def two-arg-ior logior bignum-logical-ior
)
977 (def two-arg-xor logxor bignum-logical-xor
)
978 ;; BIGNUM-LOGICAL-{AND,IOR,XOR} need not return a bignum, so must
979 ;; call the generic LOGNOT...
980 (def two-arg-eqv logeqv
(lambda (x y
) (lognot (bignum-logical-xor x y
))))
982 (lambda (x y
) (lognot (bignum-logical-and x y
)))
983 #!+sb-doc
"Complement the logical AND of INTEGER1 and INTEGER2.")
985 (lambda (x y
) (lognot (bignum-logical-ior x y
)))
986 #!+sb-doc
"Complement the logical AND of INTEGER1 and INTEGER2.")
987 ;; ... but BIGNUM-LOGICAL-NOT on a bignum will always return a bignum
988 (def logandc1 logandc1
989 (lambda (x y
) (bignum-logical-and (bignum-logical-not x
) y
))
990 #!+sb-doc
"Bitwise AND (LOGNOT INTEGER1) with INTEGER2.")
991 (def logandc2 logandc2
992 (lambda (x y
) (bignum-logical-and x
(bignum-logical-not y
)))
993 #!+sb-doc
"Bitwise AND INTEGER1 with (LOGNOT INTEGER2).")
995 (lambda (x y
) (bignum-logical-ior (bignum-logical-not x
) y
))
996 #!+sb-doc
"Bitwise OR (LOGNOT INTEGER1) with INTEGER2.")
998 (lambda (x y
) (bignum-logical-ior x
(bignum-logical-not y
)))
999 #!+sb-doc
"Bitwise OR INTEGER1 with (LOGNOT INTEGER2)."))
1001 (defun logcount (integer)
1003 "Count the number of 1 bits if INTEGER is positive, and the number of 0 bits
1004 if INTEGER is negative."
1007 (logcount (truly-the (integer 0
1008 #.
(max sb
!xc
:most-positive-fixnum
1009 (lognot sb
!xc
:most-negative-fixnum
)))
1010 (if (minusp (truly-the fixnum integer
))
1011 (lognot (truly-the fixnum integer
))
1014 (bignum-logcount integer
))))
1016 (defun logtest (integer1 integer2
)
1018 "Predicate which returns T if logand of integer1 and integer2 is not zero."
1019 (logtest integer1 integer2
))
1021 (defun logbitp (index integer
)
1023 "Predicate returns T if bit index of integer is a 1."
1024 (number-dispatch ((index integer
) (integer integer
))
1025 ((fixnum fixnum
) (if (> index
#.
(- sb
!vm
:n-word-bits sb
!vm
:n-lowtag-bits
))
1027 (not (zerop (logand integer
(ash 1 index
))))))
1028 ((fixnum bignum
) (bignum-logbitp index integer
))
1029 ((bignum (foreach fixnum bignum
)) (minusp integer
))))
1031 (defun ash (integer count
)
1033 "Shifts integer left by count places preserving sign. - count shifts right."
1034 (declare (integer integer count
))
1037 (cond ((zerop integer
)
1040 (let ((length (integer-length (truly-the fixnum integer
)))
1041 (count (truly-the fixnum count
)))
1042 (declare (fixnum length count
))
1043 (cond ((and (plusp count
)
1045 (integer-length most-positive-fixnum
)))
1046 (bignum-ashift-left (make-small-bignum integer
) count
))
1049 (ash (truly-the fixnum integer
) count
))))))
1051 (if (minusp integer
) -
1 0))
1053 (bignum-ashift-left (make-small-bignum integer
) count
))))
1056 (bignum-ashift-left integer count
)
1057 (bignum-ashift-right integer
(- count
))))))
1059 (defun integer-length (integer)
1061 "Return the number of significant bits in the absolute value of integer."
1064 (integer-length (truly-the fixnum integer
)))
1066 (bignum-integer-length integer
))))
1068 ;;;; BYTE, bytespecs, and related operations
1070 (defun byte (size position
)
1072 "Return a byte specifier which may be used by other byte functions
1074 (byte size position
))
1076 (defun byte-size (bytespec)
1078 "Return the size part of the byte specifier bytespec."
1079 (byte-size bytespec
))
1081 (defun byte-position (bytespec)
1083 "Return the position part of the byte specifier bytespec."
1084 (byte-position bytespec
))
1086 (defun ldb (bytespec integer
)
1088 "Extract the specified byte from integer, and right justify result."
1089 (ldb bytespec integer
))
1091 (defun ldb-test (bytespec integer
)
1093 "Return T if any of the specified bits in integer are 1's."
1094 (ldb-test bytespec integer
))
1096 (defun mask-field (bytespec integer
)
1098 "Extract the specified byte from integer, but do not right justify result."
1099 (mask-field bytespec integer
))
1101 (defun dpb (newbyte bytespec integer
)
1103 "Return new integer with newbyte in specified position, newbyte is right justified."
1104 (dpb newbyte bytespec integer
))
1106 (defun deposit-field (newbyte bytespec integer
)
1108 "Return new integer with newbyte in specified position, newbyte is not right justified."
1109 (deposit-field newbyte bytespec integer
))
1111 (defun %ldb
(size posn integer
)
1112 (logand (ash integer
(- posn
))
1115 (defun %mask-field
(size posn integer
)
1116 (logand integer
(ash (1- (ash 1 size
)) posn
)))
1118 (defun %dpb
(newbyte size posn integer
)
1119 (let ((mask (1- (ash 1 size
))))
1120 (logior (logand integer
(lognot (ash mask posn
)))
1121 (ash (logand newbyte mask
) posn
))))
1123 (defun %deposit-field
(newbyte size posn integer
)
1124 (let ((mask (ash (ldb (byte size
0) -
1) posn
)))
1125 (logior (logand newbyte mask
)
1126 (logand integer
(lognot mask
)))))
1130 ;;; The boole function dispaches to any logic operation depending on
1131 ;;; the value of a variable. Presently, legal selector values are [0..15].
1132 ;;; boole is open coded for calls with a constant selector. or with calls
1133 ;;; using any of the constants declared below.
1135 (defconstant boole-clr
0
1137 "Boole function op, makes BOOLE return 0.")
1139 (defconstant boole-set
1
1141 "Boole function op, makes BOOLE return -1.")
1143 (defconstant boole-1
2
1145 "Boole function op, makes BOOLE return integer1.")
1147 (defconstant boole-2
3
1149 "Boole function op, makes BOOLE return integer2.")
1151 (defconstant boole-c1
4
1153 "Boole function op, makes BOOLE return complement of integer1.")
1155 (defconstant boole-c2
5
1157 "Boole function op, makes BOOLE return complement of integer2.")
1159 (defconstant boole-and
6
1161 "Boole function op, makes BOOLE return logand of integer1 and integer2.")
1163 (defconstant boole-ior
7
1165 "Boole function op, makes BOOLE return logior of integer1 and integer2.")
1167 (defconstant boole-xor
8
1169 "Boole function op, makes BOOLE return logxor of integer1 and integer2.")
1171 (defconstant boole-eqv
9
1173 "Boole function op, makes BOOLE return logeqv of integer1 and integer2.")
1175 (defconstant boole-nand
10
1177 "Boole function op, makes BOOLE return log nand of integer1 and integer2.")
1179 (defconstant boole-nor
11
1181 "Boole function op, makes BOOLE return lognor of integer1 and integer2.")
1183 (defconstant boole-andc1
12
1185 "Boole function op, makes BOOLE return logandc1 of integer1 and integer2.")
1187 (defconstant boole-andc2
13
1189 "Boole function op, makes BOOLE return logandc2 of integer1 and integer2.")
1191 (defconstant boole-orc1
14
1193 "Boole function op, makes BOOLE return logorc1 of integer1 and integer2.")
1195 (defconstant boole-orc2
15
1197 "Boole function op, makes BOOLE return logorc2 of integer1 and integer2.")
1199 (defun boole (op integer1 integer2
)
1201 "Bit-wise boolean function on two integers. Function chosen by OP:
1219 (0 (boole 0 integer1 integer2
))
1220 (1 (boole 1 integer1 integer2
))
1221 (2 (boole 2 integer1 integer2
))
1222 (3 (boole 3 integer1 integer2
))
1223 (4 (boole 4 integer1 integer2
))
1224 (5 (boole 5 integer1 integer2
))
1225 (6 (boole 6 integer1 integer2
))
1226 (7 (boole 7 integer1 integer2
))
1227 (8 (boole 8 integer1 integer2
))
1228 (9 (boole 9 integer1 integer2
))
1229 (10 (boole 10 integer1 integer2
))
1230 (11 (boole 11 integer1 integer2
))
1231 (12 (boole 12 integer1 integer2
))
1232 (13 (boole 13 integer1 integer2
))
1233 (14 (boole 14 integer1 integer2
))
1234 (15 (boole 15 integer1 integer2
))
1235 (t (error 'type-error
:datum op
:expected-type
'(mod 16)))))
1239 (defun gcd (&rest numbers
)
1241 "Return the greatest common divisor of the arguments, which must be
1242 integers. Gcd with no arguments is defined to be 0."
1243 (cond ((null numbers
) 0)
1244 ((null (cdr numbers
)) (abs (the integer
(car numbers
))))
1246 (do ((gcd (the integer
(car numbers
))
1247 (gcd gcd
(the integer
(car rest
))))
1248 (rest (cdr numbers
) (cdr rest
)))
1250 (declare (integer gcd
)
1253 (defun lcm (&rest numbers
)
1255 "Return the least common multiple of one or more integers. LCM of no
1256 arguments is defined to be 1."
1257 (cond ((null numbers
) 1)
1258 ((null (cdr numbers
)) (abs (the integer
(car numbers
))))
1260 (do ((lcm (the integer
(car numbers
))
1261 (lcm lcm
(the integer
(car rest
))))
1262 (rest (cdr numbers
) (cdr rest
)))
1264 (declare (integer lcm
) (list rest
))))))
1266 (defun two-arg-lcm (n m
)
1267 (declare (integer n m
))
1268 (if (or (zerop n
) (zerop m
))
1270 ;; KLUDGE: I'm going to assume that it was written this way
1271 ;; originally for a reason. However, this is a somewhat
1272 ;; complicated way of writing the algorithm in the CLHS page for
1273 ;; LCM, and I don't know why. To be investigated. -- CSR,
1277 (multiple-value-bind (max min
)
1281 (* (truncate max
(gcd n m
)) min
)))))
1283 ;;; Do the GCD of two integer arguments. With fixnum arguments, we use the
1284 ;;; binary GCD algorithm from Knuth's seminumerical algorithms (slightly
1285 ;;; structurified), otherwise we call BIGNUM-GCD. We pick off the special case
1286 ;;; of 0 before the dispatch so that the bignum code doesn't have to worry
1287 ;;; about "small bignum" zeros.
1288 (defun two-arg-gcd (u v
)
1289 (cond ((eql u
0) (abs v
))
1292 (number-dispatch ((u integer
) (v integer
))
1295 (declare (optimize (speed 3) (safety 0)))
1297 (u (abs u
) (ash u -
1))
1298 (v (abs v
) (ash v -
1)))
1299 ((oddp (logior u v
))
1300 (do ((temp (if (oddp u
) (- v
) (ash u -
1))
1303 (declare (fixnum temp
))
1310 (let ((res (ash u k
)))
1311 (declare (type (signed-byte 31) res
)
1312 (optimize (inhibit-warnings 3)))
1314 (declare (type (mod 30) k
)
1315 (type (signed-byte 31) u v
)))))
1319 (bignum-gcd u
(make-small-bignum v
)))
1321 (bignum-gcd (make-small-bignum u
) v
))))))
1323 ;;; From discussion on comp.lang.lisp and Akira Kurihara.
1326 "Return the root of the nearest integer less than n which is a perfect
1328 (declare (type unsigned-byte n
) (values unsigned-byte
))
1329 ;; Theoretically (> n 7), i.e., n-len-quarter > 0.
1330 (if (and (fixnump n
) (<= n
24))
1336 (let* ((n-len-quarter (ash (integer-length n
) -
2))
1337 (n-half (ash n
(- (ash n-len-quarter
1))))
1338 (n-half-isqrt (isqrt n-half
))
1339 (init-value (ash (1+ n-half-isqrt
) n-len-quarter
)))
1341 (let ((iterated-value
1342 (ash (+ init-value
(truncate n init-value
)) -
1)))
1343 (unless (< iterated-value init-value
)
1344 (return init-value
))
1345 (setq init-value iterated-value
))))))
1347 ;;;; miscellaneous number predicates
1349 (macrolet ((def (name doc
)
1350 `(defun ,name
(number) ,doc
(,name number
))))
1351 (def zerop
"Is this number zero?")
1352 (def plusp
"Is this real number strictly positive?")
1353 (def minusp
"Is this real number strictly negative?")
1354 (def oddp
"Is this integer odd?")
1355 (def evenp
"Is this integer even?"))
1357 ;;;; modular functions
1360 (flet ((definition (name lambda-list width pattern
)
1361 (assert (sb!xc
:subtypep
`(unsigned-byte ,width
)
1362 'bignum-element-type
))
1363 `(defun ,name
,lambda-list
1364 (flet ((prepare-argument (x)
1365 (declare (integer x
))
1367 ((unsigned-byte ,width
) x
)
1368 (bignum-element-type (logand x
,pattern
))
1369 (fixnum (logand x
,pattern
))
1370 (bignum (logand (%bignum-ref x
0) ,pattern
)))))
1371 (,name
,@(loop for arg in lambda-list
1372 collect
`(prepare-argument ,arg
)))))))
1373 (loop for infos being each hash-value of sb
!c
::*modular-funs
*
1374 ;; FIXME: We need to process only "toplevel" functions
1375 unless
(eq infos
:good
)
1376 do
(loop for info in infos
1377 for name
= (sb!c
::modular-fun-info-name info
)
1378 and width
= (sb!c
::modular-fun-info-width info
)
1379 and lambda-list
= (sb!c
::modular-fun-info-lambda-list info
)
1380 for pattern
= (1- (ash 1 width
))
1381 do
(forms (definition name lambda-list width pattern
)))))
