1 ;;;; This file contains floating-point-specific transforms, and may be
2 ;;;; somewhat implementation-dependent in its assumptions of what the
5 ;;;; This software is part of the SBCL system. See the README file for
8 ;;;; This software is derived from the CMU CL system, which was
9 ;;;; written at Carnegie Mellon University and released into the
10 ;;;; public domain. The software is in the public domain and is
11 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
12 ;;;; files for more information.
18 (defknown %single-float
(real) single-float
20 (defknown %double-float
(real) double-float
23 (deftransform float
((n f
) (* single-float
) *)
26 (deftransform float
((n f
) (* double-float
) *)
29 (deftransform float
((n) *)
34 (deftransform %single-float
((n) (single-float) *)
37 (deftransform %double-float
((n) (double-float) *)
41 (macrolet ((frob (fun type
)
42 `(deftransform random
((num &optional state
)
43 (,type
&optional
*) *)
44 "Use inline float operations."
45 '(,fun num
(or state
*random-state
*)))))
46 (frob %random-single-float single-float
)
47 (frob %random-double-float double-float
))
49 ;;; Return an expression to generate an integer of N-BITS many random
50 ;;; bits, using the minimal number of random chunks possible.
51 (defun generate-random-expr-for-power-of-2 (n-bits state
)
52 (declare (type (integer 1 #.sb
!vm
:n-word-bits
) n-bits
))
53 (multiple-value-bind (n-chunk-bits chunk-expr
)
54 (cond ((<= n-bits n-random-chunk-bits
)
55 (values n-random-chunk-bits
`(random-chunk ,state
)))
56 ((<= n-bits
(* 2 n-random-chunk-bits
))
57 (values (* 2 n-random-chunk-bits
) `(big-random-chunk ,state
)))
59 (error "Unexpectedly small N-RANDOM-CHUNK-BITS")))
60 (if (< n-bits n-chunk-bits
)
61 `(logand ,(1- (ash 1 n-bits
)) ,chunk-expr
)
64 ;;; This transform for compile-time constant word-sized integers
65 ;;; generates an accept-reject loop to achieve equidistribution of the
66 ;;; returned values. Several optimizations are done: If NUM is a power
67 ;;; of two no loop is needed. If the random chunk size is half the word
68 ;;; size only one chunk is used where sufficient. For values of NUM
69 ;;; where it is possible and results in faster code, the rejection
70 ;;; probability is reduced by accepting all values below the largest
71 ;;; multiple of the limit that fits into one or two chunks and and doing
72 ;;; a division to get the random value into the desired range.
73 (deftransform random
((num &optional state
)
74 ((constant-arg (integer 1 #.
(expt 2 sb
!vm
:n-word-bits
)))
77 :policy
(and (> speed compilation-speed
)
79 "optimize to inlined RANDOM-CHUNK operations"
80 (let ((num (lvar-value num
)))
83 (flet ((chunk-n-bits-and-expr (n-bits)
84 (cond ((<= n-bits n-random-chunk-bits
)
85 (values n-random-chunk-bits
86 '(random-chunk (or state
*random-state
*))))
87 ((<= n-bits
(* 2 n-random-chunk-bits
))
88 (values (* 2 n-random-chunk-bits
)
89 '(big-random-chunk (or state
*random-state
*))))
91 (error "Unexpectedly small N-RANDOM-CHUNK-BITS")))))
92 (if (zerop (logand num
(1- num
)))
93 ;; NUM is a power of 2.
94 (let ((n-bits (integer-length (1- num
))))
95 (multiple-value-bind (n-chunk-bits chunk-expr
)
96 (chunk-n-bits-and-expr n-bits
)
97 (if (< n-bits n-chunk-bits
)
98 `(logand ,(1- (ash 1 n-bits
)) ,chunk-expr
)
100 ;; Generate an accept-reject loop.
101 (let ((n-bits (integer-length num
)))
102 (multiple-value-bind (n-chunk-bits chunk-expr
)
103 (chunk-n-bits-and-expr n-bits
)
104 (if (or (> (* num
3) (expt 2 n-chunk-bits
))
105 (logbitp (- n-bits
2) num
))
106 ;; Division can't help as the quotient is below 3,
107 ;; or is too costly as the rejection probability
108 ;; without it is already small (namely at most 1/4
109 ;; with the given test, which is experimentally a
110 ;; reasonable threshold and cheap to test for).
112 (let ((bits ,(generate-random-expr-for-power-of-2
113 n-bits
'(or state
*random-state
*))))
116 (let ((d (truncate (expt 2 n-chunk-bits
) num
)))
118 (let ((bits ,chunk-expr
))
119 (when (< bits
,(* num d
))
120 (return (values (truncate bits
,d
)))))))))))))))
125 (defknown make-single-float
((signed-byte 32)) single-float
128 (defknown make-double-float
((signed-byte 32) (unsigned-byte 32)) double-float
132 (deftransform make-single-float
((bits)
134 "Conditional constant folding"
135 (unless (constant-lvar-p bits
)
136 (give-up-ir1-transform))
137 (let* ((bits (lvar-value bits
))
138 (float (make-single-float bits
)))
139 (when (float-nan-p float
)
140 (give-up-ir1-transform))
144 (deftransform make-double-float
((hi lo
)
145 ((signed-byte 32) (unsigned-byte 32)))
146 "Conditional constant folding"
147 (unless (and (constant-lvar-p hi
)
148 (constant-lvar-p lo
))
149 (give-up-ir1-transform))
150 (let* ((hi (lvar-value hi
))
152 (float (make-double-float hi lo
)))
153 (when (float-nan-p float
)
154 (give-up-ir1-transform))
157 (defknown single-float-bits
(single-float) (signed-byte 32)
158 (movable foldable flushable
))
160 (defknown double-float-high-bits
(double-float) (signed-byte 32)
161 (movable foldable flushable
))
163 (defknown double-float-low-bits
(double-float) (unsigned-byte 32)
164 (movable foldable flushable
))
166 (deftransform float-sign
((float &optional float2
)
167 (single-float &optional single-float
) *)
169 (let ((temp (gensym)))
170 `(let ((,temp
(abs float2
)))
171 (if (minusp (single-float-bits float
)) (- ,temp
) ,temp
)))
172 '(if (minusp (single-float-bits float
)) -
1f0
1f0
)))
174 (deftransform float-sign
((float &optional float2
)
175 (double-float &optional double-float
) *)
177 (let ((temp (gensym)))
178 `(let ((,temp
(abs float2
)))
179 (if (minusp (double-float-high-bits float
)) (- ,temp
) ,temp
)))
180 '(if (minusp (double-float-high-bits float
)) -
1d0
1d0
)))
182 ;;;; DECODE-FLOAT, INTEGER-DECODE-FLOAT, and SCALE-FLOAT
184 (defknown decode-single-float
(single-float)
185 (values single-float single-float-exponent
(single-float -
1f0
1f0
))
186 (movable foldable flushable
))
188 (defknown decode-double-float
(double-float)
189 (values double-float double-float-exponent
(double-float -
1d0
1d0
))
190 (movable foldable flushable
))
192 (defknown integer-decode-single-float
(single-float)
193 (values single-float-significand single-float-int-exponent
(integer -
1 1))
194 (movable foldable flushable
))
196 (defknown integer-decode-double-float
(double-float)
197 (values double-float-significand double-float-int-exponent
(integer -
1 1))
198 (movable foldable flushable
))
200 (defknown scale-single-float
(single-float integer
) single-float
201 (movable foldable flushable
))
203 (defknown scale-double-float
(double-float integer
) double-float
204 (movable foldable flushable
))
206 (deftransform decode-float
((x) (single-float) *)
207 '(decode-single-float x
))
209 (deftransform decode-float
((x) (double-float) *)
210 '(decode-double-float x
))
212 (deftransform integer-decode-float
((x) (single-float) *)
213 '(integer-decode-single-float x
))
215 (deftransform integer-decode-float
((x) (double-float) *)
216 '(integer-decode-double-float x
))
218 (deftransform scale-float
((f ex
) (single-float *) *)
220 ((csubtypep (lvar-type ex
)
221 (specifier-type '(signed-byte 32)))
222 '(coerce (%scalbn
(coerce f
'double-float
) ex
) 'single-float
))
224 '(scale-single-float f ex
))))
226 (deftransform scale-float
((f ex
) (double-float *) *)
228 ((csubtypep (lvar-type ex
)
229 (specifier-type '(signed-byte 32)))
232 '(scale-double-float f ex
))))
234 ;;; Given a number X, create a form suitable as a bound for an
235 ;;; interval. Make the bound open if OPEN-P is T. NIL remains NIL.
236 ;;; FIXME: as this is a constructor, shouldn't it be named MAKE-BOUND?
237 #!-sb-fluid
(declaim (inline set-bound
))
238 (defun set-bound (x open-p
)
239 (if (and x open-p
) (list x
) x
))
241 ;;; What is the CROSS-FLOAT-INFINITY-KLUDGE?
243 ;;; SBCL's own implementation of floating point supports floating
244 ;;; point infinities. Some of the old CMU CL :PROPAGATE-FLOAT-TYPE and
245 ;;; :PROPAGATE-FUN-TYPE code, like the DEFOPTIMIZERs below, uses this
246 ;;; floating point support. Thus, we have to avoid running it on the
247 ;;; cross-compilation host, since we're not guaranteed that the
248 ;;; cross-compilation host will support floating point infinities.
250 ;;; If we wanted to live dangerously, we could conditionalize the code
251 ;;; with #+(OR SBCL SB-XC) instead. That way, if the cross-compilation
252 ;;; host happened to be SBCL, we'd be able to run the infinity-using
254 ;;; * SBCL itself gets built with more complete optimization.
256 ;;; * You get a different SBCL depending on what your cross-compilation
258 ;;; So far the pros and cons seem seem to be mostly academic, since
259 ;;; AFAIK (WHN 2001-08-28) the propagate-foo-type optimizations aren't
260 ;;; actually important in compiling SBCL itself. If this changes, then
261 ;;; we have to decide:
262 ;;; * Go for simplicity, leaving things as they are.
263 ;;; * Go for performance at the expense of conceptual clarity,
264 ;;; using #+(OR SBCL SB-XC) and otherwise leaving the build
266 ;;; * Go for performance at the expense of build time, using
267 ;;; #+(OR SBCL SB-XC) and also making SBCL do not just
268 ;;; make-host-1.sh and make-host-2.sh, but a third step
269 ;;; make-host-3.sh where it builds itself under itself. (Such a
270 ;;; 3-step build process could also help with other things, e.g.
271 ;;; using specialized arrays to represent debug information.)
272 ;;; * Rewrite the code so that it doesn't depend on unportable
273 ;;; floating point infinities.
275 ;;; optimizers for SCALE-FLOAT. If the float has bounds, new bounds
276 ;;; are computed for the result, if possible.
277 #-sb-xc-host
; (See CROSS-FLOAT-INFINITY-KLUDGE.)
280 (defun scale-float-derive-type-aux (f ex same-arg
)
281 (declare (ignore same-arg
))
282 (flet ((scale-bound (x n
)
283 ;; We need to be a bit careful here and catch any overflows
284 ;; that might occur. We can ignore underflows which become
288 (scale-float (type-bound-number x
) n
)
289 (floating-point-overflow ()
292 (when (and (numeric-type-p f
) (numeric-type-p ex
))
293 (let ((f-lo (numeric-type-low f
))
294 (f-hi (numeric-type-high f
))
295 (ex-lo (numeric-type-low ex
))
296 (ex-hi (numeric-type-high ex
))
300 (if (< (float-sign (type-bound-number f-hi
)) 0.0)
302 (setf new-hi
(scale-bound f-hi ex-lo
)))
304 (setf new-hi
(scale-bound f-hi ex-hi
)))))
306 (if (< (float-sign (type-bound-number f-lo
)) 0.0)
308 (setf new-lo
(scale-bound f-lo ex-hi
)))
310 (setf new-lo
(scale-bound f-lo ex-lo
)))))
311 (make-numeric-type :class
(numeric-type-class f
)
312 :format
(numeric-type-format f
)
316 (defoptimizer (scale-single-float derive-type
) ((f ex
))
317 (two-arg-derive-type f ex
#'scale-float-derive-type-aux
318 #'scale-single-float
))
319 (defoptimizer (scale-double-float derive-type
) ((f ex
))
320 (two-arg-derive-type f ex
#'scale-float-derive-type-aux
321 #'scale-double-float
))
323 ;;; DEFOPTIMIZERs for %SINGLE-FLOAT and %DOUBLE-FLOAT. This makes the
324 ;;; FLOAT function return the correct ranges if the input has some
325 ;;; defined range. Quite useful if we want to convert some type of
326 ;;; bounded integer into a float.
328 ((frob (fun type most-negative most-positive
)
329 (let ((aux-name (symbolicate fun
"-DERIVE-TYPE-AUX")))
331 (defun ,aux-name
(num)
332 ;; When converting a number to a float, the limits are
334 (let* ((lo (bound-func (lambda (x)
335 (if (< x
,most-negative
)
338 (numeric-type-low num
)
340 (hi (bound-func (lambda (x)
341 (if (< ,most-positive x
)
344 (numeric-type-high num
)
346 (specifier-type `(,',type
,(or lo
'*) ,(or hi
'*)))))
348 (defoptimizer (,fun derive-type
) ((num))
350 (one-arg-derive-type num
#',aux-name
#',fun
)
353 (frob %single-float single-float
354 most-negative-single-float most-positive-single-float
)
355 (frob %double-float double-float
356 most-negative-double-float most-positive-double-float
))
361 (defun safe-ctype-for-single-coercion-p (x)
362 ;; See comment in SAFE-SINGLE-COERCION-P -- this deals with the same
363 ;; problem, but in the context of evaluated and compiled (+ <int> <single>)
364 ;; giving different result if we fail to check for this.
365 (or (not (csubtypep x
(specifier-type 'integer
)))
367 (csubtypep x
(specifier-type `(integer ,most-negative-exactly-single-float-fixnum
368 ,most-positive-exactly-single-float-fixnum
)))
370 (csubtypep x
(specifier-type 'fixnum
))))
372 ;;; Do some stuff to recognize when the loser is doing mixed float and
373 ;;; rational arithmetic, or different float types, and fix it up. If
374 ;;; we don't, he won't even get so much as an efficiency note.
375 (deftransform float-contagion-arg1
((x y
) * * :defun-only t
:node node
)
376 (if (or (not (types-equal-or-intersect (lvar-type y
) (specifier-type 'single-float
)))
377 (safe-ctype-for-single-coercion-p (lvar-type x
)))
378 `(,(lvar-fun-name (basic-combination-fun node
))
380 (give-up-ir1-transform)))
381 (deftransform float-contagion-arg2
((x y
) * * :defun-only t
:node node
)
382 (if (or (not (types-equal-or-intersect (lvar-type x
) (specifier-type 'single-float
)))
383 (safe-ctype-for-single-coercion-p (lvar-type y
)))
384 `(,(lvar-fun-name (basic-combination-fun node
))
386 (give-up-ir1-transform)))
388 (dolist (x '(+ * / -
))
389 (%deftransform x
'(function (rational float
) *) #'float-contagion-arg1
)
390 (%deftransform x
'(function (float rational
) *) #'float-contagion-arg2
))
392 (dolist (x '(= < > + * / -
))
393 (%deftransform x
'(function (single-float double-float
) *)
394 #'float-contagion-arg1
)
395 (%deftransform x
'(function (double-float single-float
) *)
396 #'float-contagion-arg2
))
398 (macrolet ((def (type &rest args
)
399 `(deftransform * ((x y
) (,type
(constant-arg (member ,@args
))) *
401 :policy
(zerop float-accuracy
))
402 "optimize multiplication by one"
403 (let ((y (lvar-value y
)))
407 (def single-float
1.0 -
1.0)
408 (def double-float
1.0d0 -
1.0d0
))
410 ;;; Return the reciprocal of X if it can be represented exactly, NIL otherwise.
411 (defun maybe-exact-reciprocal (x)
414 (multiple-value-bind (significand exponent sign
)
415 (integer-decode-float x
)
416 ;; only powers of 2 can be inverted exactly
417 (unless (zerop (logand significand
(1- significand
)))
418 (return-from maybe-exact-reciprocal nil
))
419 (let ((expected (/ sign significand
(expt 2 exponent
)))
421 (multiple-value-bind (significand exponent sign
)
422 (integer-decode-float reciprocal
)
423 ;; Denorms can't be inverted safely.
424 (and (eql expected
(* sign significand
(expt 2 exponent
)))
426 (error () (return-from maybe-exact-reciprocal nil
)))))
428 ;;; Replace constant division by multiplication with exact reciprocal,
430 (macrolet ((def (type)
431 `(deftransform / ((x y
) (,type
(constant-arg ,type
)) *
433 "convert to multiplication by reciprocal"
434 (let ((n (lvar-value y
)))
435 (if (policy node
(zerop float-accuracy
))
437 (let ((r (maybe-exact-reciprocal n
)))
440 (give-up-ir1-transform
441 "~S does not have an exact reciprocal"
446 ;;; Optimize addition and subtraction of zero
447 (macrolet ((def (op type
&rest args
)
448 `(deftransform ,op
((x y
) (,type
(constant-arg (member ,@args
))) *
450 :policy
(zerop float-accuracy
))
452 ;; No signed zeros, thanks.
453 (def + single-float
0 0.0)
454 (def - single-float
0 0.0)
455 (def + double-float
0 0.0 0.0d0
)
456 (def - double-float
0 0.0 0.0d0
))
458 ;;; On most platforms (+ x x) is faster than (* x 2)
459 (macrolet ((def (type &rest args
)
460 `(deftransform * ((x y
) (,type
(constant-arg (member ,@args
))))
462 (def single-float
2 2.0)
463 (def double-float
2 2.0 2.0d0
))
465 ;;; Prevent ZEROP, PLUSP, and MINUSP from losing horribly. We can't in
466 ;;; general float rational args to comparison, since Common Lisp
467 ;;; semantics says we are supposed to compare as rationals, but we can
468 ;;; do it for any rational that has a precise representation as a
469 ;;; float (such as 0).
470 (macrolet ((frob (op)
471 `(deftransform ,op
((x y
) (float rational
) *)
472 "open-code FLOAT to RATIONAL comparison"
473 (unless (constant-lvar-p y
)
474 (give-up-ir1-transform
475 "The RATIONAL value isn't known at compile time."))
476 (let ((val (lvar-value y
)))
477 (unless (eql (rational (float val
)) val
)
478 (give-up-ir1-transform
479 "~S doesn't have a precise float representation."
481 `(,',op x
(float y x
)))))
486 ;;;; irrational derive-type methods
488 ;;; Derive the result to be float for argument types in the
489 ;;; appropriate domain.
490 #+sb-xc-host
; (See CROSS-FLOAT-INFINITY-KLUDGE.)
491 (dolist (stuff '((asin (real -
1.0 1.0))
492 (acos (real -
1.0 1.0))
494 (atanh (real -
1.0 1.0))
496 (destructuring-bind (name type
) stuff
497 (let ((type (specifier-type type
)))
498 (setf (fun-info-derive-type (fun-info-or-lose name
))
500 (declare (type combination call
))
501 (when (csubtypep (lvar-type
502 (first (combination-args call
)))
504 (specifier-type 'float
)))))))
506 #+sb-xc-host
; (See CROSS-FLOAT-INFINITY-KLUDGE.)
507 (defoptimizer (log derive-type
) ((x &optional y
))
508 (when (and (csubtypep (lvar-type x
)
509 (specifier-type '(real 0.0)))
511 (csubtypep (lvar-type y
)
512 (specifier-type '(real 0.0)))))
513 (specifier-type 'float
)))
515 ;;;; irrational transforms
517 (defknown (%tan %sinh %asinh %atanh %log %logb %log10 %tan-quick
)
518 (double-float) double-float
519 (movable foldable flushable
))
521 (defknown (%sin %cos %tanh %sin-quick %cos-quick
)
522 (double-float) (double-float -
1.0d0
1.0d0
)
523 (movable foldable flushable
))
525 (defknown (%asin %atan
)
527 (double-float #.
(coerce (- (/ pi
2)) 'double-float
)
528 #.
(coerce (/ pi
2) 'double-float
))
529 (movable foldable flushable
))
532 (double-float) (double-float 0.0d0
#.
(coerce pi
'double-float
))
533 (movable foldable flushable
))
536 (double-float) (double-float 1.0d0
)
537 (movable foldable flushable
))
539 (defknown (%acosh %exp %sqrt
)
540 (double-float) (double-float 0.0d0
)
541 (movable foldable flushable
))
544 (double-float) (double-float -
1d0
)
545 (movable foldable flushable
))
548 (double-float double-float
) (double-float 0d0
)
549 (movable foldable flushable
))
552 (double-float double-float
) double-float
553 (movable foldable flushable
))
556 (double-float double-float
)
557 (double-float #.
(coerce (- pi
) 'double-float
)
558 #.
(coerce pi
'double-float
))
559 (movable foldable flushable
))
562 (double-float double-float
) double-float
563 (movable foldable flushable
))
566 (double-float (signed-byte 32)) double-float
567 (movable foldable flushable
))
570 (double-float) double-float
571 (movable foldable flushable
))
573 (macrolet ((def (name prim rtype
)
575 (deftransform ,name
((x) (single-float) ,rtype
)
576 `(coerce (,',prim
(coerce x
'double-float
)) 'single-float
))
577 (deftransform ,name
((x) (double-float) ,rtype
)
581 (def sqrt %sqrt float
)
582 (def asin %asin float
)
583 (def acos %acos float
)
589 (def acosh %acosh float
)
590 (def atanh %atanh float
))
592 ;;; The argument range is limited on the x86 FP trig. functions. A
593 ;;; post-test can detect a failure (and load a suitable result), but
594 ;;; this test is avoided if possible.
595 (macrolet ((def (name prim prim-quick
)
596 (declare (ignorable prim-quick
))
598 (deftransform ,name
((x) (single-float) *)
599 #!+x86
(cond ((csubtypep (lvar-type x
)
600 (specifier-type '(single-float
601 (#.
(- (expt 2f0
63)))
603 `(coerce (,',prim-quick
(coerce x
'double-float
))
607 "unable to avoid inline argument range check~@
608 because the argument range (~S) was not within 2^63"
609 (type-specifier (lvar-type x
)))
610 `(coerce (,',prim
(coerce x
'double-float
)) 'single-float
)))
611 #!-x86
`(coerce (,',prim
(coerce x
'double-float
)) 'single-float
))
612 (deftransform ,name
((x) (double-float) *)
613 #!+x86
(cond ((csubtypep (lvar-type x
)
614 (specifier-type '(double-float
615 (#.
(- (expt 2d0
63)))
620 "unable to avoid inline argument range check~@
621 because the argument range (~S) was not within 2^63"
622 (type-specifier (lvar-type x
)))
624 #!-x86
`(,',prim x
)))))
625 (def sin %sin %sin-quick
)
626 (def cos %cos %cos-quick
)
627 (def tan %tan %tan-quick
))
629 (deftransform atan
((x y
) (single-float single-float
) *)
630 `(coerce (%atan2
(coerce x
'double-float
) (coerce y
'double-float
))
632 (deftransform atan
((x y
) (double-float double-float
) *)
635 (deftransform expt
((x y
) ((single-float 0f0
) single-float
) *)
636 `(coerce (%pow
(coerce x
'double-float
) (coerce y
'double-float
))
638 (deftransform expt
((x y
) ((double-float 0d0
) double-float
) *)
640 (deftransform expt
((x y
) ((single-float 0f0
) (signed-byte 32)) *)
641 `(coerce (%pow
(coerce x
'double-float
) (coerce y
'double-float
))
643 (deftransform expt
((x y
) ((double-float 0d0
) (signed-byte 32)) *)
644 `(%pow x
(coerce y
'double-float
)))
646 ;;; ANSI says log with base zero returns zero.
647 (deftransform log
((x y
) (float float
) float
)
648 '(if (zerop y
) y
(/ (log x
) (log y
))))
650 ;;; Handle some simple transformations.
652 (deftransform abs
((x) ((complex double-float
)) double-float
)
653 '(%hypot
(realpart x
) (imagpart x
)))
655 (deftransform abs
((x) ((complex single-float
)) single-float
)
656 '(coerce (%hypot
(coerce (realpart x
) 'double-float
)
657 (coerce (imagpart x
) 'double-float
))
660 (deftransform phase
((x) ((complex double-float
)) double-float
)
661 '(%atan2
(imagpart x
) (realpart x
)))
663 (deftransform phase
((x) ((complex single-float
)) single-float
)
664 '(coerce (%atan2
(coerce (imagpart x
) 'double-float
)
665 (coerce (realpart x
) 'double-float
))
668 (deftransform phase
((x) ((float)) float
)
669 '(if (minusp (float-sign x
))
673 ;;; The number is of type REAL.
674 (defun numeric-type-real-p (type)
675 (and (numeric-type-p type
)
676 (eq (numeric-type-complexp type
) :real
)))
678 ;;; Coerce a numeric type bound to the given type while handling
679 ;;; exclusive bounds.
680 (defun coerce-numeric-bound (bound type
)
683 (list (coerce (car bound
) type
))
684 (coerce bound type
))))
686 #-sb-xc-host
; (See CROSS-FLOAT-INFINITY-KLUDGE.)
689 ;;;; optimizers for elementary functions
691 ;;;; These optimizers compute the output range of the elementary
692 ;;;; function, based on the domain of the input.
694 ;;; Generate a specifier for a complex type specialized to the same
695 ;;; type as the argument.
696 (defun complex-float-type (arg)
697 (declare (type numeric-type arg
))
698 (let* ((format (case (numeric-type-class arg
)
699 ((integer rational
) 'single-float
)
700 (t (numeric-type-format arg
))))
701 (float-type (or format
'float
)))
702 (specifier-type `(complex ,float-type
))))
704 ;;; Compute a specifier like '(OR FLOAT (COMPLEX FLOAT)), except float
705 ;;; should be the right kind of float. Allow bounds for the float
707 (defun float-or-complex-float-type (arg &optional lo hi
)
708 (declare (type numeric-type arg
))
709 (let* ((format (case (numeric-type-class arg
)
710 ((integer rational
) 'single-float
)
711 (t (numeric-type-format arg
))))
712 (float-type (or format
'float
))
713 (lo (coerce-numeric-bound lo float-type
))
714 (hi (coerce-numeric-bound hi float-type
)))
715 (specifier-type `(or (,float-type
,(or lo
'*) ,(or hi
'*))
716 (complex ,float-type
)))))
720 (eval-when (:compile-toplevel
:execute
)
721 ;; So the problem with this hack is that it's actually broken. If
722 ;; the host does not have long floats, then setting *R-D-F-F* to
723 ;; LONG-FLOAT doesn't actually buy us anything. FIXME.
724 (setf *read-default-float-format
*
725 #!+long-float
'long-float
#!-long-float
'double-float
))
726 ;;; Test whether the numeric-type ARG is within the domain specified by
727 ;;; DOMAIN-LOW and DOMAIN-HIGH, consider negative and positive zero to
729 #-sb-xc-host
; (See CROSS-FLOAT-INFINITY-KLUDGE.)
730 (defun domain-subtypep (arg domain-low domain-high
)
731 (declare (type numeric-type arg
)
732 (type (or real null
) domain-low domain-high
))
733 (let* ((arg-lo (numeric-type-low arg
))
734 (arg-lo-val (type-bound-number arg-lo
))
735 (arg-hi (numeric-type-high arg
))
736 (arg-hi-val (type-bound-number arg-hi
)))
737 ;; Check that the ARG bounds are correctly canonicalized.
738 (when (and arg-lo
(floatp arg-lo-val
) (zerop arg-lo-val
) (consp arg-lo
)
739 (minusp (float-sign arg-lo-val
)))
740 (compiler-notify "float zero bound ~S not correctly canonicalized?" arg-lo
)
741 (setq arg-lo
0e0 arg-lo-val arg-lo
))
742 (when (and arg-hi
(zerop arg-hi-val
) (floatp arg-hi-val
) (consp arg-hi
)
743 (plusp (float-sign arg-hi-val
)))
744 (compiler-notify "float zero bound ~S not correctly canonicalized?" arg-hi
)
745 (setq arg-hi
(ecase *read-default-float-format
*
746 (double-float (load-time-value (make-unportable-float :double-float-negative-zero
)))
748 (long-float (load-time-value (make-unportable-float :long-float-negative-zero
))))
750 (flet ((fp-neg-zero-p (f) ; Is F -0.0?
751 (and (floatp f
) (zerop f
) (minusp (float-sign f
))))
752 (fp-pos-zero-p (f) ; Is F +0.0?
753 (and (floatp f
) (zerop f
) (plusp (float-sign f
)))))
754 (and (or (null domain-low
)
755 (and arg-lo
(>= arg-lo-val domain-low
)
756 (not (and (fp-pos-zero-p domain-low
)
757 (fp-neg-zero-p arg-lo
)))))
758 (or (null domain-high
)
759 (and arg-hi
(<= arg-hi-val domain-high
)
760 (not (and (fp-neg-zero-p domain-high
)
761 (fp-pos-zero-p arg-hi
)))))))))
762 (eval-when (:compile-toplevel
:execute
)
763 (setf *read-default-float-format
* 'single-float
))
765 ;;; The basic interval type. It can handle open and closed intervals.
766 ;;; A bound is open if it is a list containing a number, just like
767 ;;; Lisp says. NIL means unbounded.
768 (defstruct (interval (:constructor %make-interval
(low high
))
772 #-sb-xc-host
; (See CROSS-FLOAT-INFINITY-KLUDGE.)
775 ;;; Handle monotonic functions of a single variable whose domain is
776 ;;; possibly part of the real line. ARG is the variable, FUN is the
777 ;;; function, and DOMAIN is a specifier that gives the (real) domain
778 ;;; of the function. If ARG is a subset of the DOMAIN, we compute the
779 ;;; bounds directly. Otherwise, we compute the bounds for the
780 ;;; intersection between ARG and DOMAIN, and then append a complex
781 ;;; result, which occurs for the parts of ARG not in the DOMAIN.
783 ;;; Negative and positive zero are considered distinct within
784 ;;; DOMAIN-LOW and DOMAIN-HIGH.
786 ;;; DEFAULT-LOW and DEFAULT-HIGH are the lower and upper bounds if we
787 ;;; can't compute the bounds using FUN.
788 (defun elfun-derive-type-simple (arg fun domain-low domain-high
789 default-low default-high
790 &optional
(increasingp t
))
791 (declare (type (or null real
) domain-low domain-high
))
794 (cond ((eq (numeric-type-complexp arg
) :complex
)
795 (complex-float-type arg
))
796 ((numeric-type-real-p arg
)
797 ;; The argument is real, so let's find the intersection
798 ;; between the argument and the domain of the function.
799 ;; We compute the bounds on the intersection, and for
800 ;; everything else, we return a complex number of the
802 (multiple-value-bind (intersection difference
)
803 (interval-intersection/difference
(numeric-type->interval arg
)
809 ;; Process the intersection.
810 (let* ((low (interval-low intersection
))
811 (high (interval-high intersection
))
812 (res-lo (or (bound-func fun
(if increasingp low high
) nil
)
814 (res-hi (or (bound-func fun
(if increasingp high low
) nil
)
816 (format (case (numeric-type-class arg
)
817 ((integer rational
) 'single-float
)
818 (t (numeric-type-format arg
))))
819 (bound-type (or format
'float
))
824 :low
(coerce-numeric-bound res-lo bound-type
)
825 :high
(coerce-numeric-bound res-hi bound-type
))))
826 ;; If the ARG is a subset of the domain, we don't
827 ;; have to worry about the difference, because that
829 (if (or (null difference
)
830 ;; Check whether the arg is within the domain.
831 (domain-subtypep arg domain-low domain-high
))
834 (specifier-type `(complex ,bound-type
))))))
836 ;; No intersection so the result must be purely complex.
837 (complex-float-type arg
)))))
839 (float-or-complex-float-type arg default-low default-high
))))))
842 ((frob (name domain-low domain-high def-low-bnd def-high-bnd
843 &key
(increasingp t
))
844 (let ((num (gensym)))
845 `(defoptimizer (,name derive-type
) ((,num
))
849 (elfun-derive-type-simple arg
#',name
850 ,domain-low
,domain-high
851 ,def-low-bnd
,def-high-bnd
854 ;; These functions are easy because they are defined for the whole
856 (frob exp nil nil
0 nil
)
857 (frob sinh nil nil nil nil
)
858 (frob tanh nil nil -
1 1)
859 (frob asinh nil nil nil nil
)
861 ;; These functions are only defined for part of the real line. The
862 ;; condition selects the desired part of the line.
863 (frob asin -
1d0
1d0
(- (/ pi
2)) (/ pi
2))
864 ;; Acos is monotonic decreasing, so we need to swap the function
865 ;; values at the lower and upper bounds of the input domain.
866 (frob acos -
1d0
1d0
0 pi
:increasingp nil
)
867 (frob acosh
1d0 nil nil nil
)
868 (frob atanh -
1d0
1d0 -
1 1)
869 ;; Kahan says that (sqrt -0.0) is -0.0, so use a specifier that
871 (frob sqrt
(load-time-value (make-unportable-float :double-float-negative-zero
)) nil
0 nil
))
873 ;;; Compute bounds for (expt x y). This should be easy since (expt x
874 ;;; y) = (exp (* y (log x))). However, computations done this way
875 ;;; have too much roundoff. Thus we have to do it the hard way.
876 (defun safe-expt (x y
)
878 (when (< (abs y
) 10000)
883 ;;; Handle the case when x >= 1.
884 (defun interval-expt-> (x y
)
885 (case (interval-range-info y
0d0
)
887 ;; Y is positive and log X >= 0. The range of exp(y * log(x)) is
888 ;; obviously non-negative. We just have to be careful for
889 ;; infinite bounds (given by nil).
890 (let ((lo (safe-expt (type-bound-number (interval-low x
))
891 (type-bound-number (interval-low y
))))
892 (hi (safe-expt (type-bound-number (interval-high x
))
893 (type-bound-number (interval-high y
)))))
894 (list (make-interval :low
(or lo
1) :high hi
))))
896 ;; Y is negative and log x >= 0. The range of exp(y * log(x)) is
897 ;; obviously [0, 1]. However, underflow (nil) means 0 is the
899 (let ((lo (safe-expt (type-bound-number (interval-high x
))
900 (type-bound-number (interval-low y
))))
901 (hi (safe-expt (type-bound-number (interval-low x
))
902 (type-bound-number (interval-high y
)))))
903 (list (make-interval :low
(or lo
0) :high
(or hi
1)))))
905 ;; Split the interval in half.
906 (destructuring-bind (y- y
+)
907 (interval-split 0 y t
)
908 (list (interval-expt-> x y-
)
909 (interval-expt-> x y
+))))))
911 ;;; Handle the case when x <= 1
912 (defun interval-expt-< (x y
)
913 (case (interval-range-info x
0d0
)
915 ;; The case of 0 <= x <= 1 is easy
916 (case (interval-range-info y
)
918 ;; Y is positive and log X <= 0. The range of exp(y * log(x)) is
919 ;; obviously [0, 1]. We just have to be careful for infinite bounds
921 (let ((lo (safe-expt (type-bound-number (interval-low x
))
922 (type-bound-number (interval-high y
))))
923 (hi (safe-expt (type-bound-number (interval-high x
))
924 (type-bound-number (interval-low y
)))))
925 (list (make-interval :low
(or lo
0) :high
(or hi
1)))))
927 ;; Y is negative and log x <= 0. The range of exp(y * log(x)) is
928 ;; obviously [1, inf].
929 (let ((hi (safe-expt (type-bound-number (interval-low x
))
930 (type-bound-number (interval-low y
))))
931 (lo (safe-expt (type-bound-number (interval-high x
))
932 (type-bound-number (interval-high y
)))))
933 (list (make-interval :low
(or lo
1) :high hi
))))
935 ;; Split the interval in half
936 (destructuring-bind (y- y
+)
937 (interval-split 0 y t
)
938 (list (interval-expt-< x y-
)
939 (interval-expt-< x y
+))))))
941 ;; The case where x <= 0. Y MUST be an INTEGER for this to work!
942 ;; The calling function must insure this! For now we'll just
943 ;; return the appropriate unbounded float type.
944 (list (make-interval :low nil
:high nil
)))
946 (destructuring-bind (neg pos
)
947 (interval-split 0 x t t
)
948 (list (interval-expt-< neg y
)
949 (interval-expt-< pos y
))))))
951 ;;; Compute bounds for (expt x y).
952 (defun interval-expt (x y
)
953 (case (interval-range-info x
1)
956 (interval-expt-> x y
))
959 (interval-expt-< x y
))
961 (destructuring-bind (left right
)
962 (interval-split 1 x t t
)
963 (list (interval-expt left y
)
964 (interval-expt right y
))))))
966 (defun fixup-interval-expt (bnd x-int y-int x-type y-type
)
967 (declare (ignore x-int
))
968 ;; Figure out what the return type should be, given the argument
969 ;; types and bounds and the result type and bounds.
970 (cond ((csubtypep x-type
(specifier-type 'integer
))
971 ;; an integer to some power
972 (case (numeric-type-class y-type
)
974 ;; Positive integer to an integer power is either an
975 ;; integer or a rational.
976 (let ((lo (or (interval-low bnd
) '*))
977 (hi (or (interval-high bnd
) '*)))
978 (if (and (interval-low y-int
)
979 (>= (type-bound-number (interval-low y-int
)) 0))
980 (specifier-type `(integer ,lo
,hi
))
981 (specifier-type `(rational ,lo
,hi
)))))
983 ;; Positive integer to rational power is either a rational
984 ;; or a single-float.
985 (let* ((lo (interval-low bnd
))
986 (hi (interval-high bnd
))
988 (floor (type-bound-number lo
))
991 (ceiling (type-bound-number hi
))
993 (f-lo (or (bound-func #'float lo nil
)
995 (f-hi (or (bound-func #'float hi nil
)
997 (specifier-type `(or (rational ,int-lo
,int-hi
)
998 (single-float ,f-lo
, f-hi
)))))
1000 ;; A positive integer to a float power is a float.
1001 (let ((format (numeric-type-format y-type
)))
1003 (modified-numeric-type
1005 :low
(coerce-numeric-bound (interval-low bnd
) format
)
1006 :high
(coerce-numeric-bound (interval-high bnd
) format
))))
1008 ;; A positive integer to a number is a number (for now).
1009 (specifier-type 'number
))))
1010 ((csubtypep x-type
(specifier-type 'rational
))
1011 ;; a rational to some power
1012 (case (numeric-type-class y-type
)
1014 ;; A positive rational to an integer power is always a rational.
1015 (specifier-type `(rational ,(or (interval-low bnd
) '*)
1016 ,(or (interval-high bnd
) '*))))
1018 ;; A positive rational to rational power is either a rational
1019 ;; or a single-float.
1020 (let* ((lo (interval-low bnd
))
1021 (hi (interval-high bnd
))
1023 (floor (type-bound-number lo
))
1026 (ceiling (type-bound-number hi
))
1028 (f-lo (or (bound-func #'float lo nil
)
1030 (f-hi (or (bound-func #'float hi nil
)
1032 (specifier-type `(or (rational ,int-lo
,int-hi
)
1033 (single-float ,f-lo
, f-hi
)))))
1035 ;; A positive rational to a float power is a float.
1036 (let ((format (numeric-type-format y-type
)))
1038 (modified-numeric-type
1040 :low
(coerce-numeric-bound (interval-low bnd
) format
)
1041 :high
(coerce-numeric-bound (interval-high bnd
) format
))))
1043 ;; A positive rational to a number is a number (for now).
1044 (specifier-type 'number
))))
1045 ((csubtypep x-type
(specifier-type 'float
))
1046 ;; a float to some power
1047 (case (numeric-type-class y-type
)
1048 ((or integer rational
)
1049 ;; A positive float to an integer or rational power is
1051 (let ((format (numeric-type-format x-type
)))
1056 :low
(coerce-numeric-bound (interval-low bnd
) format
)
1057 :high
(coerce-numeric-bound (interval-high bnd
) format
))))
1059 ;; A positive float to a float power is a float of the
1061 (let ((format (float-format-max (numeric-type-format x-type
)
1062 (numeric-type-format y-type
))))
1067 :low
(coerce-numeric-bound (interval-low bnd
) format
)
1068 :high
(coerce-numeric-bound (interval-high bnd
) format
))))
1070 ;; A positive float to a number is a number (for now)
1071 (specifier-type 'number
))))
1073 ;; A number to some power is a number.
1074 (specifier-type 'number
))))
1076 (defun merged-interval-expt (x y
)
1077 (let* ((x-int (numeric-type->interval x
))
1078 (y-int (numeric-type->interval y
)))
1079 (mapcar (lambda (type)
1080 (fixup-interval-expt type x-int y-int x y
))
1081 (flatten-list (interval-expt x-int y-int
)))))
1083 (defun expt-derive-type-aux (x y same-arg
)
1084 (declare (ignore same-arg
))
1085 (cond ((or (not (numeric-type-real-p x
))
1086 (not (numeric-type-real-p y
)))
1087 ;; Use numeric contagion if either is not real.
1088 (numeric-contagion x y
))
1089 ((csubtypep y
(specifier-type 'integer
))
1090 ;; A real raised to an integer power is well-defined.
1091 (merged-interval-expt x y
))
1092 ;; A real raised to a non-integral power can be a float or a
1094 ((or (csubtypep x
(specifier-type '(rational 0)))
1095 (csubtypep x
(specifier-type '(float (0d0)))))
1096 ;; But a positive real to any power is well-defined.
1097 (merged-interval-expt x y
))
1098 ((and (csubtypep x
(specifier-type 'rational
))
1099 (csubtypep y
(specifier-type 'rational
)))
1100 ;; A rational to the power of a rational could be a rational
1101 ;; or a possibly-complex single float
1102 (specifier-type '(or rational single-float
(complex single-float
))))
1104 ;; a real to some power. The result could be a real or a
1106 (float-or-complex-float-type (numeric-contagion x y
)))))
1108 (defoptimizer (expt derive-type
) ((x y
))
1109 (two-arg-derive-type x y
#'expt-derive-type-aux
#'expt
))
1111 ;;; Note we must assume that a type including 0.0 may also include
1112 ;;; -0.0 and thus the result may be complex -infinity + i*pi.
1113 (defun log-derive-type-aux-1 (x)
1114 (elfun-derive-type-simple x
#'log
0d0 nil nil nil
))
1116 (defun log-derive-type-aux-2 (x y same-arg
)
1117 (let ((log-x (log-derive-type-aux-1 x
))
1118 (log-y (log-derive-type-aux-1 y
))
1119 (accumulated-list nil
))
1120 ;; LOG-X or LOG-Y might be union types. We need to run through
1121 ;; the union types ourselves because /-DERIVE-TYPE-AUX doesn't.
1122 (dolist (x-type (prepare-arg-for-derive-type log-x
))
1123 (dolist (y-type (prepare-arg-for-derive-type log-y
))
1124 (push (/-derive-type-aux x-type y-type same-arg
) accumulated-list
)))
1125 (apply #'type-union
(flatten-list accumulated-list
))))
1127 (defoptimizer (log derive-type
) ((x &optional y
))
1129 (two-arg-derive-type x y
#'log-derive-type-aux-2
#'log
)
1130 (one-arg-derive-type x
#'log-derive-type-aux-1
#'log
)))
1132 (defun atan-derive-type-aux-1 (y)
1133 (elfun-derive-type-simple y
#'atan nil nil
(- (/ pi
2)) (/ pi
2)))
1135 (defun atan-derive-type-aux-2 (y x same-arg
)
1136 (declare (ignore same-arg
))
1137 ;; The hard case with two args. We just return the max bounds.
1138 (let ((result-type (numeric-contagion y x
)))
1139 (cond ((and (numeric-type-real-p x
)
1140 (numeric-type-real-p y
))
1141 (let* (;; FIXME: This expression for FORMAT seems to
1142 ;; appear multiple times, and should be factored out.
1143 (format (case (numeric-type-class result-type
)
1144 ((integer rational
) 'single-float
)
1145 (t (numeric-type-format result-type
))))
1146 (bound-format (or format
'float
)))
1147 (make-numeric-type :class
'float
1150 :low
(coerce (- pi
) bound-format
)
1151 :high
(coerce pi bound-format
))))
1153 ;; The result is a float or a complex number
1154 (float-or-complex-float-type result-type
)))))
1156 (defoptimizer (atan derive-type
) ((y &optional x
))
1158 (two-arg-derive-type y x
#'atan-derive-type-aux-2
#'atan
)
1159 (one-arg-derive-type y
#'atan-derive-type-aux-1
#'atan
)))
1161 (defun cosh-derive-type-aux (x)
1162 ;; We note that cosh x = cosh |x| for all real x.
1163 (elfun-derive-type-simple
1164 (if (numeric-type-real-p x
)
1165 (abs-derive-type-aux x
)
1167 #'cosh nil nil
0 nil
))
1169 (defoptimizer (cosh derive-type
) ((num))
1170 (one-arg-derive-type num
#'cosh-derive-type-aux
#'cosh
))
1172 (defun phase-derive-type-aux (arg)
1173 (let* ((format (case (numeric-type-class arg
)
1174 ((integer rational
) 'single-float
)
1175 (t (numeric-type-format arg
))))
1176 (bound-type (or format
'float
)))
1177 (cond ((numeric-type-real-p arg
)
1178 (case (interval-range-info (numeric-type->interval arg
) 0.0)
1180 ;; The number is positive, so the phase is 0.
1181 (make-numeric-type :class
'float
1184 :low
(coerce 0 bound-type
)
1185 :high
(coerce 0 bound-type
)))
1187 ;; The number is always negative, so the phase is pi.
1188 (make-numeric-type :class
'float
1191 :low
(coerce pi bound-type
)
1192 :high
(coerce pi bound-type
)))
1194 ;; We can't tell. The result is 0 or pi. Use a union
1197 (make-numeric-type :class
'float
1200 :low
(coerce 0 bound-type
)
1201 :high
(coerce 0 bound-type
))
1202 (make-numeric-type :class
'float
1205 :low
(coerce pi bound-type
)
1206 :high
(coerce pi bound-type
))))))
1208 ;; We have a complex number. The answer is the range -pi
1209 ;; to pi. (-pi is included because we have -0.)
1210 (make-numeric-type :class
'float
1213 :low
(coerce (- pi
) bound-type
)
1214 :high
(coerce pi bound-type
))))))
1216 (defoptimizer (phase derive-type
) ((num))
1217 (one-arg-derive-type num
#'phase-derive-type-aux
#'phase
))
1221 (deftransform realpart
((x) ((complex rational
)) *)
1223 (deftransform imagpart
((x) ((complex rational
)) *)
1226 ;;; Make REALPART and IMAGPART return the appropriate types. This
1227 ;;; should help a lot in optimized code.
1228 (defun realpart-derive-type-aux (type)
1229 (let ((class (numeric-type-class type
))
1230 (format (numeric-type-format type
)))
1231 (cond ((numeric-type-real-p type
)
1232 ;; The realpart of a real has the same type and range as
1234 (make-numeric-type :class class
1237 :low
(numeric-type-low type
)
1238 :high
(numeric-type-high type
)))
1240 ;; We have a complex number. The result has the same type
1241 ;; as the real part, except that it's real, not complex,
1243 (make-numeric-type :class class
1246 :low
(numeric-type-low type
)
1247 :high
(numeric-type-high type
))))))
1249 (defoptimizer (realpart derive-type
) ((num))
1250 (one-arg-derive-type num
#'realpart-derive-type-aux
#'realpart
))
1252 (defun imagpart-derive-type-aux (type)
1253 (let ((class (numeric-type-class type
))
1254 (format (numeric-type-format type
)))
1255 (cond ((numeric-type-real-p type
)
1256 ;; The imagpart of a real has the same type as the input,
1257 ;; except that it's zero.
1258 (let ((bound-format (or format class
'real
)))
1259 (make-numeric-type :class class
1262 :low
(coerce 0 bound-format
)
1263 :high
(coerce 0 bound-format
))))
1265 ;; We have a complex number. The result has the same type as
1266 ;; the imaginary part, except that it's real, not complex,
1268 (make-numeric-type :class class
1271 :low
(numeric-type-low type
)
1272 :high
(numeric-type-high type
))))))
1274 (defoptimizer (imagpart derive-type
) ((num))
1275 (one-arg-derive-type num
#'imagpart-derive-type-aux
#'imagpart
))
1277 (defun complex-derive-type-aux-1 (re-type)
1278 (if (numeric-type-p re-type
)
1279 (make-numeric-type :class
(numeric-type-class re-type
)
1280 :format
(numeric-type-format re-type
)
1281 :complexp
(if (csubtypep re-type
1282 (specifier-type 'rational
))
1285 :low
(numeric-type-low re-type
)
1286 :high
(numeric-type-high re-type
))
1287 (specifier-type 'complex
)))
1289 (defun complex-derive-type-aux-2 (re-type im-type same-arg
)
1290 (declare (ignore same-arg
))
1291 (if (and (numeric-type-p re-type
)
1292 (numeric-type-p im-type
))
1293 ;; Need to check to make sure numeric-contagion returns the
1294 ;; right type for what we want here.
1296 ;; Also, what about rational canonicalization, like (complex 5 0)
1297 ;; is 5? So, if the result must be complex, we make it so.
1298 ;; If the result might be complex, which happens only if the
1299 ;; arguments are rational, we make it a union type of (or
1300 ;; rational (complex rational)).
1301 (let* ((element-type (numeric-contagion re-type im-type
))
1302 (rat-result-p (csubtypep element-type
1303 (specifier-type 'rational
))))
1305 (type-union element-type
1307 `(complex ,(numeric-type-class element-type
))))
1308 (make-numeric-type :class
(numeric-type-class element-type
)
1309 :format
(numeric-type-format element-type
)
1310 :complexp
(if rat-result-p
1313 (specifier-type 'complex
)))
1315 #-sb-xc-host
; (See CROSS-FLOAT-INFINITY-KLUDGE.)
1316 (defoptimizer (complex derive-type
) ((re &optional im
))
1318 (two-arg-derive-type re im
#'complex-derive-type-aux-2
#'complex
)
1319 (one-arg-derive-type re
#'complex-derive-type-aux-1
#'complex
)))
1321 ;;; Define some transforms for complex operations. We do this in lieu
1322 ;;; of complex operation VOPs.
1323 (macrolet ((frob (type)
1325 (deftransform complex
((r) (,type
))
1326 '(complex r
,(coerce 0 type
)))
1327 (deftransform complex
((r i
) (,type
(and real
(not ,type
))))
1328 '(complex r
(truly-the ,type
(coerce i
',type
))))
1329 (deftransform complex
((r i
) ((and real
(not ,type
)) ,type
))
1330 '(complex (truly-the ,type
(coerce r
',type
)) i
))
1332 #!-complex-float-vops
1333 (deftransform %negate
((z) ((complex ,type
)) *)
1334 '(complex (%negate
(realpart z
)) (%negate
(imagpart z
))))
1335 ;; complex addition and subtraction
1336 #!-complex-float-vops
1337 (deftransform + ((w z
) ((complex ,type
) (complex ,type
)) *)
1338 '(complex (+ (realpart w
) (realpart z
))
1339 (+ (imagpart w
) (imagpart z
))))
1340 #!-complex-float-vops
1341 (deftransform -
((w z
) ((complex ,type
) (complex ,type
)) *)
1342 '(complex (- (realpart w
) (realpart z
))
1343 (- (imagpart w
) (imagpart z
))))
1344 ;; Add and subtract a complex and a real.
1345 #!-complex-float-vops
1346 (deftransform + ((w z
) ((complex ,type
) real
) *)
1347 `(complex (+ (realpart w
) z
)
1348 (+ (imagpart w
) ,(coerce 0 ',type
))))
1349 #!-complex-float-vops
1350 (deftransform + ((z w
) (real (complex ,type
)) *)
1351 `(complex (+ (realpart w
) z
)
1352 (+ (imagpart w
) ,(coerce 0 ',type
))))
1353 ;; Add and subtract a real and a complex number.
1354 #!-complex-float-vops
1355 (deftransform -
((w z
) ((complex ,type
) real
) *)
1356 `(complex (- (realpart w
) z
)
1357 (- (imagpart w
) ,(coerce 0 ',type
))))
1358 #!-complex-float-vops
1359 (deftransform -
((z w
) (real (complex ,type
)) *)
1360 `(complex (- z
(realpart w
))
1361 (- ,(coerce 0 ',type
) (imagpart w
))))
1362 ;; Multiply and divide two complex numbers.
1363 #!-complex-float-vops
1364 (deftransform * ((x y
) ((complex ,type
) (complex ,type
)) *)
1365 '(let* ((rx (realpart x
))
1369 (complex (- (* rx ry
) (* ix iy
))
1370 (+ (* rx iy
) (* ix ry
)))))
1371 (deftransform / ((x y
) ((complex ,type
) (complex ,type
)) *)
1372 #!-complex-float-vops
1373 '(let* ((rx (realpart x
))
1377 (if (> (abs ry
) (abs iy
))
1378 (let* ((r (/ iy ry
))
1379 (dn (+ ry
(* r iy
))))
1380 (complex (/ (+ rx
(* ix r
)) dn
)
1381 (/ (- ix
(* rx r
)) dn
)))
1382 (let* ((r (/ ry iy
))
1383 (dn (+ iy
(* r ry
))))
1384 (complex (/ (+ (* rx r
) ix
) dn
)
1385 (/ (- (* ix r
) rx
) dn
)))))
1386 #!+complex-float-vops
1387 `(let* ((cs (conjugate (sb!vm
::swap-complex x
)))
1390 (if (> (abs ry
) (abs iy
))
1391 (let* ((r (/ iy ry
))
1392 (dn (+ ry
(* r iy
))))
1393 (/ (+ x
(* cs r
)) dn
))
1394 (let* ((r (/ ry iy
))
1395 (dn (+ iy
(* r ry
))))
1396 (/ (+ (* x r
) cs
) dn
)))))
1397 ;; Multiply a complex by a real or vice versa.
1398 #!-complex-float-vops
1399 (deftransform * ((w z
) ((complex ,type
) real
) *)
1400 '(complex (* (realpart w
) z
) (* (imagpart w
) z
)))
1401 #!-complex-float-vops
1402 (deftransform * ((z w
) (real (complex ,type
)) *)
1403 '(complex (* (realpart w
) z
) (* (imagpart w
) z
)))
1404 ;; Divide a complex by a real or vice versa.
1405 #!-complex-float-vops
1406 (deftransform / ((w z
) ((complex ,type
) real
) *)
1407 '(complex (/ (realpart w
) z
) (/ (imagpart w
) z
)))
1408 (deftransform / ((x y
) (,type
(complex ,type
)) *)
1409 #!-complex-float-vops
1410 '(let* ((ry (realpart y
))
1412 (if (> (abs ry
) (abs iy
))
1413 (let* ((r (/ iy ry
))
1414 (dn (+ ry
(* r iy
))))
1416 (/ (- (* x r
)) dn
)))
1417 (let* ((r (/ ry iy
))
1418 (dn (+ iy
(* r ry
))))
1419 (complex (/ (* x r
) dn
)
1421 #!+complex-float-vops
1422 '(let* ((ry (realpart y
))
1424 (if (> (abs ry
) (abs iy
))
1425 (let* ((r (/ iy ry
))
1426 (dn (+ ry
(* r iy
))))
1427 (/ (complex x
(- (* x r
))) dn
))
1428 (let* ((r (/ ry iy
))
1429 (dn (+ iy
(* r ry
))))
1430 (/ (complex (* x r
) (- x
)) dn
)))))
1431 ;; conjugate of complex number
1432 #!-complex-float-vops
1433 (deftransform conjugate
((z) ((complex ,type
)) *)
1434 '(complex (realpart z
) (- (imagpart z
))))
1436 (deftransform cis
((z) ((,type
)) *)
1437 '(complex (cos z
) (sin z
)))
1439 #!-complex-float-vops
1440 (deftransform = ((w z
) ((complex ,type
) (complex ,type
)) *)
1441 '(and (= (realpart w
) (realpart z
))
1442 (= (imagpart w
) (imagpart z
))))
1443 #!-complex-float-vops
1444 (deftransform = ((w z
) ((complex ,type
) real
) *)
1445 '(and (= (realpart w
) z
) (zerop (imagpart w
))))
1446 #!-complex-float-vops
1447 (deftransform = ((w z
) (real (complex ,type
)) *)
1448 '(and (= (realpart z
) w
) (zerop (imagpart z
)))))))
1451 (frob double-float
))
1453 ;;; Here are simple optimizers for SIN, COS, and TAN. They do not
1454 ;;; produce a minimal range for the result; the result is the widest
1455 ;;; possible answer. This gets around the problem of doing range
1456 ;;; reduction correctly but still provides useful results when the
1457 ;;; inputs are union types.
1458 #-sb-xc-host
; (See CROSS-FLOAT-INFINITY-KLUDGE.)
1460 (defun trig-derive-type-aux (arg domain fun
1461 &optional def-lo def-hi
(increasingp t
))
1464 (cond ((eq (numeric-type-complexp arg
) :complex
)
1465 (make-numeric-type :class
(numeric-type-class arg
)
1466 :format
(numeric-type-format arg
)
1467 :complexp
:complex
))
1468 ((numeric-type-real-p arg
)
1469 (let* ((format (case (numeric-type-class arg
)
1470 ((integer rational
) 'single-float
)
1471 (t (numeric-type-format arg
))))
1472 (bound-type (or format
'float
)))
1473 ;; If the argument is a subset of the "principal" domain
1474 ;; of the function, we can compute the bounds because
1475 ;; the function is monotonic. We can't do this in
1476 ;; general for these periodic functions because we can't
1477 ;; (and don't want to) do the argument reduction in
1478 ;; exactly the same way as the functions themselves do
1480 (if (csubtypep arg domain
)
1481 (let ((res-lo (bound-func fun
(numeric-type-low arg
) nil
))
1482 (res-hi (bound-func fun
(numeric-type-high arg
) nil
)))
1484 (rotatef res-lo res-hi
))
1488 :low
(coerce-numeric-bound res-lo bound-type
)
1489 :high
(coerce-numeric-bound res-hi bound-type
)))
1493 :low
(and def-lo
(coerce def-lo bound-type
))
1494 :high
(and def-hi
(coerce def-hi bound-type
))))))
1496 (float-or-complex-float-type arg def-lo def-hi
))))))
1498 (defoptimizer (sin derive-type
) ((num))
1499 (one-arg-derive-type
1502 ;; Derive the bounds if the arg is in [-pi/2, pi/2].
1503 (trig-derive-type-aux
1505 (specifier-type `(float ,(- (/ pi
2)) ,(/ pi
2)))
1510 (defoptimizer (cos derive-type
) ((num))
1511 (one-arg-derive-type
1514 ;; Derive the bounds if the arg is in [0, pi].
1515 (trig-derive-type-aux arg
1516 (specifier-type `(float 0d0
,pi
))
1522 (defoptimizer (tan derive-type
) ((num))
1523 (one-arg-derive-type
1526 ;; Derive the bounds if the arg is in [-pi/2, pi/2].
1527 (trig-derive-type-aux arg
1528 (specifier-type `(float ,(- (/ pi
2)) ,(/ pi
2)))
1533 (defoptimizer (conjugate derive-type
) ((num))
1534 (one-arg-derive-type num
1536 (flet ((most-negative-bound (l h
)
1538 (if (< (type-bound-number l
) (- (type-bound-number h
)))
1540 (set-bound (- (type-bound-number h
)) (consp h
)))))
1541 (most-positive-bound (l h
)
1543 (if (> (type-bound-number h
) (- (type-bound-number l
)))
1545 (set-bound (- (type-bound-number l
)) (consp l
))))))
1546 (if (numeric-type-real-p arg
)
1548 (let ((low (numeric-type-low arg
))
1549 (high (numeric-type-high arg
)))
1550 (let ((new-low (most-negative-bound low high
))
1551 (new-high (most-positive-bound low high
)))
1552 (modified-numeric-type arg
:low new-low
:high new-high
))))))
1555 (defoptimizer (cis derive-type
) ((num))
1556 (one-arg-derive-type num
1559 `(complex ,(or (numeric-type-format arg
) 'float
))))
1564 ;;;; TRUNCATE, FLOOR, CEILING, and ROUND
1566 (macrolet ((define-frobs (fun ufun
)
1568 (defknown ,ufun
(real) integer
(movable foldable flushable
))
1569 (deftransform ,fun
((x &optional by
)
1571 (constant-arg (member 1))))
1572 '(let ((res (,ufun x
)))
1573 (values res
(- x res
)))))))
1574 (define-frobs truncate %unary-truncate
)
1575 (define-frobs round %unary-round
))
1577 (deftransform %unary-truncate
((x) (single-float))
1578 `(%unary-truncate
/single-float x
))
1579 (deftransform %unary-truncate
((x) (double-float))
1580 `(%unary-truncate
/double-float x
))
1582 ;;; Convert (TRUNCATE x y) to the obvious implementation.
1584 ;;; ...plus hair: Insert explicit coercions to appropriate float types: Python
1585 ;;; is reluctant it generate explicit integer->float coercions due to
1586 ;;; precision issues (see SAFE-SINGLE-COERCION-P &co), but this is not an
1587 ;;; issue here as there is no DERIVE-TYPE optimizer on specialized versions of
1588 ;;; %UNARY-TRUNCATE, so the derived type of TRUNCATE remains the best we can
1589 ;;; do here -- which is fine. Also take care not to add unnecassary division
1590 ;;; or multiplication by 1, since we are not able to always eliminate them,
1591 ;;; depending on FLOAT-ACCURACY. Finally, leave out the secondary value when
1592 ;;; we know it is unused: COERCE is not flushable.
1593 (macrolet ((def (type other-float-arg-types
)
1594 (let ((unary (symbolicate "%UNARY-TRUNCATE/" type
))
1595 (coerce (symbolicate "%" type
)))
1596 `(deftransform truncate
((x &optional y
)
1598 &optional
(or ,type
,@other-float-arg-types integer
))
1600 (let* ((result-type (and result
1601 (lvar-derived-type result
)))
1602 (compute-all (and (values-type-p result-type
)
1603 (not (type-single-value-p result-type
)))))
1605 (and (constant-lvar-p y
) (= 1 (lvar-value y
))))
1607 `(let ((res (,',unary x
)))
1608 (values res
(- x
(,',coerce res
))))
1609 `(let ((res (,',unary x
)))
1610 ;; Dummy secondary value!
1613 `(let* ((f (,',coerce y
))
1614 (res (,',unary
(/ x f
))))
1615 (values res
(- x
(* f
(,',coerce res
)))))
1616 `(let* ((f (,',coerce y
))
1617 (res (,',unary
(/ x f
))))
1618 ;; Dummy secondary value!
1619 (values res x
)))))))))
1620 (def single-float
())
1621 (def double-float
(single-float)))
1623 (defknown %unary-ftruncate
(real) float
(movable foldable flushable
))
1624 (defknown %unary-ftruncate
/single
(single-float) single-float
1625 (movable foldable flushable
))
1626 (defknown %unary-ftruncate
/double
(double-float) double-float
1627 (movable foldable flushable
))
1630 (defun %unary-ftruncate
/single
(x)
1631 (declare (muffle-conditions t
))
1632 (declare (type single-float x
))
1633 (declare (optimize speed
(safety 0)))
1634 (let* ((bits (single-float-bits x
))
1635 (exp (ldb sb
!vm
:single-float-exponent-byte bits
))
1636 (biased (the single-float-exponent
1637 (- exp sb
!vm
:single-float-bias
))))
1638 (declare (type (signed-byte 32) bits
))
1640 ((= exp sb
!vm
:single-float-normal-exponent-max
) x
)
1641 ((<= biased
0) (* x
0f0
))
1642 ((>= biased
(float-digits x
)) x
)
1644 (let ((frac-bits (- (float-digits x
) biased
)))
1645 (setf bits
(logandc2 bits
(- (ash 1 frac-bits
) 1)))
1646 (make-single-float bits
))))))
1649 (defun %unary-ftruncate
/double
(x)
1650 (declare (muffle-conditions t
))
1651 (declare (type double-float x
))
1652 (declare (optimize speed
(safety 0)))
1653 (let* ((high (double-float-high-bits x
))
1654 (low (double-float-low-bits x
))
1655 (exp (ldb sb
!vm
:double-float-exponent-byte high
))
1656 (biased (the double-float-exponent
1657 (- exp sb
!vm
:double-float-bias
))))
1658 (declare (type (signed-byte 32) high
)
1659 (type (unsigned-byte 32) low
))
1661 ((= exp sb
!vm
:double-float-normal-exponent-max
) x
)
1662 ((<= biased
0) (* x
0d0
))
1663 ((>= biased
(float-digits x
)) x
)
1665 (let ((frac-bits (- (float-digits x
) biased
)))
1666 (cond ((< frac-bits
32)
1667 (setf low
(logandc2 low
(- (ash 1 frac-bits
) 1))))
1670 (setf high
(logandc2 high
(- (ash 1 (- frac-bits
32)) 1)))))
1671 (make-double-float high low
))))))
1674 ((def (float-type fun
)
1675 `(deftransform %unary-ftruncate
((x) (,float-type
))
1677 (def single-float %unary-ftruncate
/single
)
1678 (def double-float %unary-ftruncate
/double
))