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 (deftransform float
((n f
) (t single-float
) *)
21 (deftransform float
((n f
) (t double-float
) *)
24 (deftransform float
((n) *)
29 (deftransform %single-float
((n) (single-float) * :important nil
)
32 (deftransform %double-float
((n) (double-float) * :important nil
)
35 (deftransform %single-float
((n) (ratio) * :important nil
)
36 '(sb-kernel::single-float-ratio n
))
38 (deftransform %double-float
((n) (ratio) * :important nil
)
39 '(sb-kernel::double-float-ratio n
))
41 (macrolet ((def (type from-type
)
42 `(deftransform ,(symbolicate "%" type
) ((n) ((or ,type
,from-type
)) * :important nil
)
43 (when (or (csubtypep (lvar-type n
) (specifier-type ',type
))
44 (csubtypep (lvar-type n
) (specifier-type ',from-type
)))
45 (give-up-ir1-transform))
46 `(if (,',(symbolicate type
"-P") n
)
48 (,',(symbolicate "%" type
) (truly-the ,',from-type n
))))))
49 (def single-float double-float
)
50 (def single-float sb-vm
:signed-word
)
51 (def single-float word
)
52 (def double-float single-float
)
53 (def double-float sb-vm
:signed-word
)
54 (def double-float word
))
56 (macrolet ((def (type from-type inline-type
)
57 `(deftransform ,(symbolicate "%" type
) ((n) (,from-type
) * :important nil
)
58 (when (or (csubtypep (lvar-type n
) (specifier-type ',inline-type
))
59 (not (types-equal-or-intersect (lvar-type n
) (specifier-type ',inline-type
))))
60 (give-up-ir1-transform))
61 '(if (typep n
',inline-type
)
62 (,(symbolicate "%" type
) (truly-the ,inline-type n
))
63 (,(symbolicate "%" type
) (truly-the (not ,inline-type
) n
))))))
64 (def single-float integer fixnum
)
65 (def double-float integer fixnum
))
68 (macrolet ((frob (fun type
)
69 `(deftransform random
((num &optional state
)
70 (,type
&optional t
) *)
71 "Use inline float operations."
72 '(,fun num
(or state
*random-state
*)))))
73 (frob %random-single-float single-float
)
74 (frob %random-double-float double-float
))
76 ;;; Return an expression to generate an integer of N-BITS many random
77 ;;; bits, using the minimal number of random chunks possible.
78 (defun generate-random-expr-for-power-of-2 (n-bits state
)
79 (declare (type (integer 1 #.sb-vm
:n-word-bits
) n-bits
))
80 (multiple-value-bind (n-chunk-bits chunk-expr
)
81 (cond ((<= n-bits n-random-chunk-bits
)
82 (values n-random-chunk-bits
`(random-chunk ,state
)))
83 ((<= n-bits
(* 2 n-random-chunk-bits
))
84 (values (* 2 n-random-chunk-bits
) `(big-random-chunk ,state
)))
86 (error "Unexpectedly small N-RANDOM-CHUNK-BITS")))
87 (if (< n-bits n-chunk-bits
)
88 `(logand ,(1- (ash 1 n-bits
)) ,chunk-expr
)
91 ;;; This transform for compile-time constant word-sized integers
92 ;;; generates an accept-reject loop to achieve equidistribution of the
93 ;;; returned values. Several optimizations are done: If NUM is a power
94 ;;; of two no loop is needed. If the random chunk size is half the word
95 ;;; size only one chunk is used where sufficient. For values of NUM
96 ;;; where it is possible and results in faster code, the rejection
97 ;;; probability is reduced by accepting all values below the largest
98 ;;; multiple of the limit that fits into one or two chunks and and doing
99 ;;; a division to get the random value into the desired range.
100 (deftransform random
((num &optional state
)
101 ((constant-arg (integer 1 #.
(expt 2 sb-vm
:n-word-bits
)))
104 :policy
(and (> speed compilation-speed
)
106 "optimize to inlined RANDOM-CHUNK operations"
107 (let ((num (lvar-value num
)))
110 (flet ((chunk-n-bits-and-expr (n-bits)
111 (cond ((<= n-bits n-random-chunk-bits
)
112 (values n-random-chunk-bits
113 '(random-chunk (or state
*random-state
*))))
114 ((<= n-bits
(* 2 n-random-chunk-bits
))
115 (values (* 2 n-random-chunk-bits
)
116 '(big-random-chunk (or state
*random-state
*))))
118 (error "Unexpectedly small N-RANDOM-CHUNK-BITS")))))
119 (if (zerop (logand num
(1- num
)))
120 ;; NUM is a power of 2.
121 (let ((n-bits (integer-length (1- num
))))
122 (multiple-value-bind (n-chunk-bits chunk-expr
)
123 (chunk-n-bits-and-expr n-bits
)
124 (if (< n-bits n-chunk-bits
)
125 `(logand ,(1- (ash 1 n-bits
)) ,chunk-expr
)
127 ;; Generate an accept-reject loop.
128 (let ((n-bits (integer-length num
)))
129 (multiple-value-bind (n-chunk-bits chunk-expr
)
130 (chunk-n-bits-and-expr n-bits
)
131 (if (or (> (* num
3) (expt 2 n-chunk-bits
))
132 (logbitp (- n-bits
2) num
))
133 ;; Division can't help as the quotient is below 3,
134 ;; or is too costly as the rejection probability
135 ;; without it is already small (namely at most 1/4
136 ;; with the given test, which is experimentally a
137 ;; reasonable threshold and cheap to test for).
139 (let ((bits ,(generate-random-expr-for-power-of-2
140 n-bits
'(or state
*random-state
*))))
143 (let ((d (truncate (expt 2 n-chunk-bits
) num
)))
145 (let ((bits ,chunk-expr
))
146 (when (< bits
,(* num d
))
147 (return (values (truncate bits
,d
)))))))))))))))
152 ;;; NaNs can not be constructed from constant bits mainly due to compiler problems
153 ;;; in so doing. See https://bugs.launchpad.net/sbcl/+bug/486812
154 (deftransform make-single-float
((bits) ((constant-arg t
)))
155 "Conditional constant folding"
156 (let ((float (make-single-float (lvar-value bits
))))
157 (if (float-nan-p float
) (give-up-ir1-transform) float
)))
159 (deftransform make-double-float
((hi lo
) ((constant-arg t
) (constant-arg t
)))
160 "Conditional constant folding"
161 (let ((float (make-double-float (lvar-value hi
) (lvar-value lo
))))
162 (if (float-nan-p float
) (give-up-ir1-transform) float
)))
164 ;;; I'd like to transition all the 64-bit backends to use the single-arg
165 ;;; %MAKE-DOUBLE-FLOAT constructor instead of the 2-arg MAKE-DOUBLE-FLOAT.
166 ;;; So we need a transform to fold constant calls for either.
168 (deftransform %make-double-float
((bits) ((constant-arg t
)))
169 "Conditional constant folding"
170 (let ((float (%make-double-float
(lvar-value bits
))))
171 (if (float-nan-p float
) (give-up-ir1-transform) float
)))
173 ;;; On the face of it, these transforms are ridiculous because if we're going
174 ;;; to express (MINUSP X) as (MINUSP (foo-FLOAT-BITS X)), then why not _always_
175 ;;; transform MINUSP of a float into an integer comparison instead of a
176 ;;; floating-point comparison, and then express this as (if (minusp float) ...)
177 ;;; rather than (if (minusp (bits float)) ...) ?
178 ;;; I suspect that the difference is that FLOAT-SIGN must remain silent
179 ;;; when given a signaling NaN.
180 (deftransform float-sign
((float &optional float2
)
181 (single-float &optional single-float
) *)
182 (if (vop-existsp :translate single-float-copysign
)
184 `(single-float-copysign float float2
)
185 `(single-float-sign float
))
187 (let ((temp (gensym)))
188 `(let ((,temp
(abs float2
)))
189 (if (minusp (single-float-bits float
)) (- ,temp
) ,temp
)))
190 '(if (minusp (single-float-bits float
)) -
1f0
1f0
))))
192 (deftransform float-sign
((float &optional float2
)
193 (double-float &optional double-float
) *)
194 ;; If words are 64 bits, then it's actually simpler to extract _all_ bits
195 ;; instead of only the upper bits.
196 (let ((bits #+64-bit
'(double-float-bits float
)
197 #-
64-bit
'(double-float-high-bits float
)))
199 (let ((temp (gensym)))
200 `(let ((,temp
(abs float2
)))
201 (if (minusp ,bits
) (- ,temp
) ,temp
)))
202 `(if (minusp ,bits
) -
1d0
1d0
))))
204 (deftransform float-sign-bit
((x) (single-float) *)
205 `(logand (ash (single-float-bits x
) -
31) 1))
206 (deftransform float-sign-bit
((x) (double-float) *)
207 #-
64-bit
`(logand (ash (double-float-high-bits x
) -
31) 1)
208 #+64-bit
`(ash (logand (double-float-bits x
) most-positive-word
) -
63))
210 (deftransform float-sign-bit-set-p
((x) (single-float) *)
211 `(logbitp 31 (single-float-bits x
)))
212 (deftransform float-sign-bit-set-p
((x) (double-float) *)
213 #-
64-bit
`(logbitp 31 (double-float-high-bits x
))
214 #+64-bit
`(logbitp 63 (double-float-bits x
)))
216 ;;; This doesn't deal with complex at the moment.
217 (deftransform signum
((x) (number))
218 (let* ((ctype (lvar-type x
))
220 (dolist (x '(single-float double-float rational
))
221 (when (csubtypep ctype
(specifier-type x
))
223 ;; SB-XC:COERCE doesn't like RATIONAL for some reason.
224 (when (eq result-type
'rational
) (setq result-type
'integer
))
227 ((plusp x
) ,(sb-xc:coerce
1 result-type
))
228 (t ,(sb-xc:coerce -
1 result-type
)))
229 (give-up-ir1-transform))))
231 ;;;; DECODE-FLOAT, INTEGER-DECODE-FLOAT, and SCALE-FLOAT
233 (defknown decode-single-float
(single-float)
234 (values single-float single-float-exponent
(member -
1f0
1f0
))
235 (movable foldable flushable
))
237 (defknown decode-double-float
(double-float)
238 (values double-float double-float-exponent
(member -
1d0
1d0
))
239 (movable foldable flushable
))
241 (defknown integer-decode-single-float
(single-float)
242 (values single-float-significand single-float-int-exponent
(member -
1 1))
243 (movable foldable flushable
))
245 (defknown integer-decode-double-float
(double-float)
246 (values double-float-significand double-float-int-exponent
(member -
1 1))
247 (movable foldable flushable
))
249 (defknown scale-single-float
(single-float integer
) single-float
250 (movable foldable flushable fixed-args unboxed-return
))
251 (defknown scale-double-float
(double-float integer
) double-float
252 (movable foldable flushable fixed-args unboxed-return
))
254 (defknown sb-kernel
::scale-single-float-maybe-overflow
255 (single-float integer
) single-float
256 (movable foldable flushable fixed-args unboxed-return
))
257 (defknown sb-kernel
::scale-single-float-maybe-underflow
258 (single-float integer
) single-float
259 (movable foldable flushable fixed-args unboxed-return
))
260 (defknown sb-kernel
::scale-double-float-maybe-overflow
261 (double-float integer
) double-float
262 (movable foldable flushable fixed-args unboxed-return
))
263 (defknown sb-kernel
::scale-double-float-maybe-underflow
264 (double-float integer
) double-float
265 (movable foldable flushable fixed-args unboxed-return
))
267 (deftransform decode-float
((x) (single-float) *)
268 '(decode-single-float x
))
270 (deftransform decode-float
((x) (double-float) *)
271 '(decode-double-float x
))
273 (deftransform integer-decode-float
((x) (single-float) *)
274 '(integer-decode-single-float x
))
276 (deftransform integer-decode-float
((x) (double-float) *)
277 '(integer-decode-double-float x
))
279 (deftransform scale-float
((f ex
) (single-float t
) *)
280 (cond #+(and x86
()) ;; this producess different results based on whether it's inlined or not
281 ((csubtypep (lvar-type ex
)
282 (specifier-type '(signed-byte 32)))
283 '(coerce (%scalbn
(coerce f
'double-float
) ex
) 'single-float
))
285 '(scale-single-float f ex
))))
287 (deftransform scale-float
((f ex
) (double-float t
) *)
289 ((csubtypep (lvar-type ex
)
290 (specifier-type '(signed-byte 32)))
293 '(scale-double-float f ex
))))
295 ;;; Given a number X, create a form suitable as a bound for an
296 ;;; interval. Make the bound open if OPEN-P is T. NIL remains NIL.
297 ;;; FIXME: as this is a constructor, shouldn't it be named MAKE-BOUND?
298 (declaim (inline set-bound
))
299 (defun set-bound (x open-p
)
300 (if (and x open-p
) (list x
) x
))
302 ;;; optimizers for SCALE-FLOAT. If the float has bounds, new bounds
303 ;;; are computed for the result, if possible.
304 (defun scale-float-derive-type-aux (f ex same-arg
)
305 (declare (ignore same-arg
))
306 (flet ((scale-bound (x n
)
307 ;; We need to be a bit careful here and catch any overflows
308 ;; that might occur. We can ignore underflows which become
312 (scale-float (type-bound-number x
) n
)
313 (floating-point-overflow ()
316 (when (and (numeric-type-p f
) (numeric-type-p ex
))
317 (let ((f-lo (numeric-type-low f
))
318 (f-hi (numeric-type-high f
))
319 (ex-lo (numeric-type-low ex
))
320 (ex-hi (numeric-type-high ex
))
324 (if (sb-xc:< (float-sign (type-bound-number f-hi
)) 0.0)
326 (setf new-hi
(scale-bound f-hi ex-lo
)))
328 (setf new-hi
(scale-bound f-hi ex-hi
)))))
330 (if (sb-xc:< (float-sign (type-bound-number f-lo
)) 0.0)
332 (setf new-lo
(scale-bound f-lo ex-hi
)))
334 (setf new-lo
(scale-bound f-lo ex-lo
)))))
335 (make-numeric-type :class
(numeric-type-class f
)
336 :format
(numeric-type-format f
)
340 (defoptimizer (scale-single-float derive-type
) ((f ex
))
341 (two-arg-derive-type f ex
#'scale-float-derive-type-aux
342 #'scale-single-float
))
343 (defoptimizer (scale-double-float derive-type
) ((f ex
))
344 (two-arg-derive-type f ex
#'scale-float-derive-type-aux
345 #'scale-double-float
))
347 ;;; DEFOPTIMIZERs for %SINGLE-FLOAT and %DOUBLE-FLOAT. This makes the
348 ;;; FLOAT function return the correct ranges if the input has some
349 ;;; defined range. Quite useful if we want to convert some type of
350 ;;; bounded integer into a float.
352 ((frob (fun type most-negative most-positive
)
353 (let ((aux-name (symbolicate fun
"-DERIVE-TYPE-AUX")))
355 (defun ,aux-name
(num)
356 ;; When converting a number to a float, the limits are
358 (let* ((lo (bound-func (lambda (x)
359 (if (sb-xc:< x
,most-negative
)
362 (numeric-type-low num
)
364 (hi (bound-func (lambda (x)
365 (if (sb-xc:< ,most-positive x
)
368 (numeric-type-high num
)
370 (specifier-type `(,',type
,(or lo
'*) ,(or hi
'*)))))
372 (defoptimizer (,fun derive-type
) ((num))
374 (one-arg-derive-type num
#',aux-name
#',fun
)
377 (frob %single-float single-float
378 most-negative-single-float most-positive-single-float
)
379 (frob %double-float double-float
380 most-negative-double-float most-positive-double-float
))
382 (defoptimizer (float derive-type
) ((number prototype
))
383 (let ((type (lvar-type prototype
)))
384 (unless (or (csubtypep type
(specifier-type 'double-float
))
385 (csubtypep type
(specifier-type 'single-float
)))
388 (one-arg-derive-type number
#'%single-float-derive-type-aux
#'%single-float
)
389 (one-arg-derive-type number
#'%double-float-derive-type-aux
#'%double-float
))
393 (macrolet ((def (type &rest args
)
394 `(deftransform * ((x y
) (,type
(constant-arg (member ,@args
))) *
396 :policy
(zerop float-accuracy
))
397 "optimize multiplication by one"
398 (let ((y (lvar-value y
)))
402 (def single-float
1.0 -
1.0)
403 (def double-float
1.0d0 -
1.0d0
))
405 ;;; Return the reciprocal of X if it can be represented exactly, NIL otherwise.
406 (defun maybe-exact-reciprocal (x)
409 (multiple-value-bind (significand exponent sign
)
410 (integer-decode-float x
)
411 ;; only powers of 2 can be inverted exactly
412 (unless (zerop (logand significand
(1- significand
)))
413 (return-from maybe-exact-reciprocal nil
))
414 (let ((expected (/ sign significand
(expt 2 exponent
)))
416 (multiple-value-bind (significand exponent sign
)
417 (integer-decode-float reciprocal
)
418 ;; Denorms can't be inverted safely.
419 (and (eql expected
(* sign significand
(expt 2 exponent
)))
421 (error () (return-from maybe-exact-reciprocal nil
)))))
423 ;;; Replace constant division by multiplication with exact reciprocal,
425 (macrolet ((def (type)
426 `(deftransform / ((x y
) (,type
(constant-arg ,type
)) *
428 "convert to multiplication by reciprocal"
429 (let ((n (lvar-value y
)))
430 (if (policy node
(zerop float-accuracy
))
432 (let ((r (maybe-exact-reciprocal n
)))
435 (give-up-ir1-transform))))))))
439 ;;; Optimize addition and subtraction of zero
440 (macrolet ((def (op type
&rest args
)
441 `(deftransform ,op
((x y
) (,type
(constant-arg (member ,@args
))) *
443 :policy
(zerop float-accuracy
))
445 ;; No signed zeros, thanks.
446 (def + single-float
0 0.0)
447 (def - single-float
0 0.0)
448 (def + double-float
0 0.0 0.0d0
)
449 (def - double-float
0 0.0 0.0d0
))
451 ;;; On most platforms (+ x x) is faster than (* x 2)
452 (macrolet ((def (type &rest args
)
453 `(deftransform * ((x y
) (,type
(constant-arg (member ,@args
))))
455 (def single-float
2 2.0)
456 (def double-float
2 2.0 2.0d0
))
458 ;;; Prevent ZEROP, PLUSP, and MINUSP from losing horribly. We can't in
459 ;;; general float rational args to comparison, since Common Lisp
460 ;;; semantics says we are supposed to compare as rationals, but we can
461 ;;; do it for any rational that has a precise representation as a
462 ;;; float (such as 0).
463 (macrolet ((frob (op &optional complex
)
464 `(deftransform ,op
((x y
) (:or
((,(if complex
465 '(complex single-float
)
469 '(complex double-float
)
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 (multiple-value-bind (low high type
)
478 (if (csubtypep (lvar-type x
) (specifier-type 'double-float
))
479 (values most-negative-double-float most-positive-double-float
481 (values most-negative-single-float most-positive-single-float
483 (unless (and (sb-xc:<= low val high
)
484 (eql (rational (coerce val type
)) val
))
485 (give-up-ir1-transform
486 "~S doesn't have a precise float representation."
488 `(,',op x
(float y
,',(if complex
499 ;;;; irrational transforms
501 (macrolet ((def (name prim rtype
)
503 (deftransform ,name
((x) (single-float) ,rtype
:node node
)
504 (delay-ir1-transform node
:ir1-phases
)
505 `(%single-float
(,',prim
(%double-float x
))))
506 (deftransform ,name
((x) (double-float) ,rtype
:node node
)
507 (delay-ir1-transform node
:ir1-phases
)
511 (def sqrt %sqrt float
)
512 (def asin %asin float
)
513 (def acos %acos float
)
519 (def acosh %acosh float
)
520 (def atanh %atanh float
))
522 ;;; The argument range is limited on the x86 FP trig. functions. A
523 ;;; post-test can detect a failure (and load a suitable result), but
524 ;;; this test is avoided if possible.
525 (macrolet ((def (name prim prim-quick
)
526 (declare (ignorable prim-quick
))
528 (deftransform ,name
((x) (single-float) *)
529 #+x86
(cond ((csubtypep (lvar-type x
)
531 `(single-float (,(sb-xc:-
(expt 2f0
63)))
533 `(coerce (,',prim-quick
(coerce x
'double-float
))
537 "unable to avoid inline argument range check~@
538 because the argument range (~S) was not within 2^63"
539 (type-specifier (lvar-type x
)))
540 `(coerce (,',prim
(coerce x
'double-float
)) 'single-float
)))
541 #-x86
`(coerce (,',prim
(coerce x
'double-float
)) 'single-float
))
542 (deftransform ,name
((x) (double-float) *)
543 #+x86
(cond ((csubtypep (lvar-type x
)
545 `(double-float (,(sb-xc:-
(expt 2d0
63)))
550 "unable to avoid inline argument range check~@
551 because the argument range (~S) was not within 2^63"
552 (type-specifier (lvar-type x
)))
554 #-x86
`(,',prim x
)))))
555 (def sin %sin %sin-quick
)
556 (def cos %cos %cos-quick
)
557 (def tan %tan %tan-quick
))
559 (deftransform atan
((x y
) (single-float single-float
) *)
560 `(coerce (%atan2
(coerce x
'double-float
) (coerce y
'double-float
))
562 (deftransform atan
((x y
) (double-float double-float
) *)
565 (deftransform expt
((x y
) (single-float single-float
) single-float
)
566 `(coerce (%pow
(coerce x
'double-float
) (coerce y
'double-float
))
568 (deftransform expt
((x y
) (double-float double-float
) double-float
)
570 (deftransform expt
((x y
) (single-float integer
) single-float
)
571 `(coerce (%pow
(coerce x
'double-float
) (coerce y
'double-float
))
573 (deftransform expt
((x y
) (double-float integer
) double-float
)
574 `(%pow x
(coerce y
'double-float
)))
576 ;;; ANSI says log with base zero returns zero.
577 (deftransform log
((x y
) (single-float single-float
) single-float
:node node
)
578 (delay-ir1-transform node
:ir1-phases
)
581 (coerce (/ (%log
(coerce x
'double-float
)) (%log
(coerce y
'double-float
)))
583 (deftransform log
((x y
) (single-float double-float
) double-float
:node node
)
584 (delay-ir1-transform node
:ir1-phases
)
587 (/ (%log
(coerce x
'double-float
)) (%log y
))))
588 (deftransform log
((x y
) (double-float single-float
) double-float
:node node
)
589 (delay-ir1-transform node
:ir1-phases
)
592 (/ (%log x
) (%log
(coerce y
'double-float
)))))
593 (deftransform log
((x y
) (double-float double-float
) double-float
:node node
)
594 (delay-ir1-transform node
:ir1-phases
)
597 (/ (%log x
) (%log y
))))
599 ;;; Handle some simple transformations.
601 (deftransform abs
((x) ((complex double-float
)) double-float
)
602 '(%hypot
(realpart x
) (imagpart x
)))
604 (deftransform abs
((x) ((complex single-float
)) single-float
)
605 '(coerce (%hypot
(coerce (realpart x
) 'double-float
)
606 (coerce (imagpart x
) 'double-float
))
609 (deftransform phase
((x) ((complex double-float
)) double-float
)
610 '(%atan2
(imagpart x
) (realpart x
)))
612 (deftransform phase
((x) ((complex single-float
)) single-float
)
613 '(coerce (%atan2
(coerce (imagpart x
) 'double-float
)
614 (coerce (realpart x
) 'double-float
))
617 (deftransform phase
((x) ((float)) float
)
618 '(if (minusp (float-sign x
))
622 ;;; The number is of type REAL.
623 (defun numeric-type-real-p (type)
624 (and (numeric-type-p type
)
625 (eq (numeric-type-complexp type
) :real
)))
627 ;;;; optimizers for elementary functions
629 ;;;; These optimizers compute the output range of the elementary
630 ;;;; function, based on the domain of the input.
632 ;;; Generate a specifier for a complex type specialized to the same
633 ;;; type as the argument.
634 (defun complex-float-type (arg)
635 (declare (type numeric-type arg
))
636 (let* ((format (case (numeric-type-class arg
)
637 ((integer rational
) 'single-float
)
638 (t (numeric-type-format arg
))))
639 (float-type (or format
'float
)))
640 (specifier-type `(complex ,float-type
))))
642 ;;; Compute a specifier like '(OR FLOAT (COMPLEX FLOAT)), except float
643 ;;; should be the right kind of float. Allow bounds for the float
645 (defun float-or-complex-float-type (arg &optional lo hi
)
647 ((numeric-type-p arg
)
648 (let* ((format (case (numeric-type-class arg
)
649 ((integer rational
) 'single-float
)
650 (t (numeric-type-format arg
))))
651 (float-type (or format
'float
))
652 (lo (coerce-numeric-bound lo float-type
))
653 (hi (coerce-numeric-bound hi float-type
)))
654 (specifier-type `(or (,float-type
,(or lo
'*) ,(or hi
'*))
655 (complex ,float-type
)))))
658 (loop for type in
(union-type-types arg
)
659 collect
(float-or-complex-float-type type lo hi
))))
660 (t (specifier-type 'number
))))
662 ;;; Test whether the numeric-type ARG is within the domain specified by
663 ;;; DOMAIN-LOW and DOMAIN-HIGH, consider negative and positive zero to
665 (defun domain-subtypep (arg domain-low domain-high
)
666 (declare (type numeric-type arg
)
667 (type (or real null
) domain-low domain-high
))
668 (let* ((arg-lo (numeric-type-low arg
))
669 (arg-lo-val (type-bound-number arg-lo
))
670 (arg-hi (numeric-type-high arg
))
671 (arg-hi-val (type-bound-number arg-hi
)))
672 ;; Check that the ARG bounds are correctly canonicalized.
673 (when (and arg-lo
(floatp arg-lo-val
) (zerop arg-lo-val
) (consp arg-lo
)
674 (minusp (float-sign arg-lo-val
)))
682 (when (and arg-hi
(zerop arg-hi-val
) (floatp arg-hi-val
) (consp arg-hi
)
683 (plusp (float-sign arg-hi-val
)))
691 (flet ((fp-neg-zero-p (f) ; Is F -0.0?
692 (and (floatp f
) (zerop f
) (float-sign-bit-set-p f
)))
693 (fp-pos-zero-p (f) ; Is F +0.0?
694 (and (floatp f
) (zerop f
) (not (float-sign-bit-set-p f
)))))
695 (and (or (null domain-low
)
696 (and arg-lo
(sb-xc:>= arg-lo-val domain-low
)
697 (not (and (fp-pos-zero-p domain-low
)
698 (fp-neg-zero-p arg-lo
)))))
699 (or (null domain-high
)
700 (and arg-hi
(sb-xc:<= arg-hi-val domain-high
)
701 (not (and (fp-neg-zero-p domain-high
)
702 (fp-pos-zero-p arg-hi
)))))))))
705 ;;; Handle monotonic functions of a single variable whose domain is
706 ;;; possibly part of the real line. ARG is the variable, FUN is the
707 ;;; function, and DOMAIN is a specifier that gives the (real) domain
708 ;;; of the function. If ARG is a subset of the DOMAIN, we compute the
709 ;;; bounds directly. Otherwise, we compute the bounds for the
710 ;;; intersection between ARG and DOMAIN, and then append a complex
711 ;;; result, which occurs for the parts of ARG not in the DOMAIN.
713 ;;; Negative and positive zero are considered distinct within
714 ;;; DOMAIN-LOW and DOMAIN-HIGH.
716 ;;; DEFAULT-LOW and DEFAULT-HIGH are the lower and upper bounds if we
717 ;;; can't compute the bounds using FUN.
718 (defun elfun-derive-type-simple (arg fun domain-low domain-high
719 default-low default-high
720 &optional
(increasingp t
))
721 (declare (type (or null real
) domain-low domain-high
))
724 (cond ((eq (numeric-type-complexp arg
) :complex
)
725 (complex-float-type arg
))
726 ((numeric-type-real-p arg
)
727 ;; The argument is real, so let's find the intersection
728 ;; between the argument and the domain of the function.
729 ;; We compute the bounds on the intersection, and for
730 ;; everything else, we return a complex number of the
732 (multiple-value-bind (intersection difference
)
733 (interval-intersection/difference
(numeric-type->interval arg
)
739 ;; Process the intersection.
740 (let* ((low (interval-low intersection
))
741 (high (interval-high intersection
))
742 (res-lo (or (bound-func fun
(if increasingp low high
) nil
)
744 (res-hi (or (bound-func fun
(if increasingp high low
) nil
)
746 (format (case (numeric-type-class arg
)
747 ((integer rational
) 'single-float
)
748 (t (numeric-type-format arg
))))
749 (bound-type (or format
'float
))
754 :low
(coerce-numeric-bound res-lo bound-type
)
755 :high
(coerce-numeric-bound res-hi bound-type
))))
756 ;; If the ARG is a subset of the domain, we don't
757 ;; have to worry about the difference, because that
759 (if (or (null difference
)
760 ;; Check whether the arg is within the domain.
761 (domain-subtypep arg domain-low domain-high
))
764 (specifier-type `(complex ,bound-type
))))))
766 ;; No intersection so the result must be purely complex.
767 (complex-float-type arg
)))))
769 (float-or-complex-float-type arg default-low default-high
))))))
772 ((frob (name domain-low domain-high def-low-bnd def-high-bnd
773 &key
(increasingp t
))
774 (let ((num (gensym)))
775 `(defoptimizer (,name derive-type
) ((,num
))
779 (elfun-derive-type-simple arg
#',name
780 ,domain-low
,domain-high
781 ,def-low-bnd
,def-high-bnd
784 ;; These functions are easy because they are defined for the whole
786 (frob exp nil nil
0 nil
)
787 (frob sinh nil nil nil nil
)
788 (frob tanh nil nil -
1 1)
789 (frob asinh nil nil nil nil
)
791 ;; These functions are only defined for part of the real line. The
792 ;; condition selects the desired part of the line.
793 (frob asin -
1d0
1d0
(sb-xc:-
(sb-xc:/ pi
2)) (sb-xc:/ pi
2))
794 ;; Acos is monotonic decreasing, so we need to swap the function
795 ;; values at the lower and upper bounds of the input domain.
796 (frob acos -
1d0
1d0
0 pi
:increasingp nil
)
797 (frob acosh
1d0 nil nil nil
)
798 (frob atanh -
1d0
1d0 -
1 1)
799 ;; Kahan says that (sqrt -0.0) is -0.0, so use a specifier that
801 (frob sqrt -
0.0d0 nil
0 nil
))
803 ;;; Compute bounds for (expt x y). This should be easy since (expt x
804 ;;; y) = (exp (* y (log x))). However, computations done this way
805 ;;; have too much roundoff. Thus we have to do it the hard way.
806 (defun safe-expt (x y
)
807 (when (and (numberp x
) (numberp y
))
809 (when (sb-xc:< (abs y
) 10000)
811 ;; Currently we can hide unanticipated errors (such as failure to use SB-XC: math
812 ;; when cross-compiling) as well as the anticipated potential problem of overflow.
813 ;; So don't handle anything when cross-compiling.
814 ;; FIXME: I think this should not handle ERROR, but just FLOATING-POINT-OVERFLOW.
816 #-sb-xc-host error
()
819 ;;; Handle the case when x >= 1.
820 (defun interval-expt-> (x y
)
821 (case (interval-range-info y
0d0
)
823 ;; Y is positive and log X >= 0. The range of exp(y * log(x)) is
824 ;; obviously non-negative. We just have to be careful for
825 ;; infinite bounds (given by nil).
826 (let ((lo (safe-expt (type-bound-number (interval-low x
))
827 (type-bound-number (interval-low y
))))
828 (hi (safe-expt (type-bound-number (interval-high x
))
829 (type-bound-number (interval-high y
)))))
830 (list (make-interval :low
(or lo
1) :high hi
))))
832 ;; Y is negative and log x >= 0. The range of exp(y * log(x)) is
833 ;; obviously [0, 1]. However, underflow (nil) means 0 is the
835 (let ((lo (safe-expt (type-bound-number (interval-high x
))
836 (type-bound-number (interval-low y
))))
837 (hi (safe-expt (type-bound-number (interval-low x
))
838 (type-bound-number (interval-high y
)))))
839 (list (make-interval :low
(or lo
0) :high
(or hi
1)))))
841 ;; Split the interval in half.
842 (destructuring-bind (y- y
+)
843 (interval-split 0 y t
)
844 (list (interval-expt-> x y-
)
845 (interval-expt-> x y
+))))))
847 ;;; Handle the case when 0 <= x <= 1
848 (defun interval-expt-< (x y
)
849 (case (interval-range-info x
0d0
)
851 ;; The case of 0 <= x <= 1 is easy
852 (case (interval-range-info y
)
854 ;; Y is positive and log X <= 0. The range of exp(y * log(x)) is
855 ;; obviously [0, 1]. We just have to be careful for infinite bounds
857 (let ((lo (safe-expt (type-bound-number (interval-low x
))
858 (type-bound-number (interval-high y
))))
859 (hi (safe-expt (type-bound-number (interval-high x
))
860 (type-bound-number (interval-low y
)))))
861 (list (make-interval :low
(or lo
0) :high
(or hi
1)))))
863 ;; Y is negative and log x <= 0. The range of exp(y * log(x)) is
864 ;; obviously [1, inf].
865 (let ((hi (safe-expt (type-bound-number (interval-low x
))
866 (type-bound-number (interval-low y
))))
867 (lo (safe-expt (type-bound-number (interval-high x
))
868 (type-bound-number (interval-high y
)))))
869 (list (make-interval :low
(or lo
1) :high hi
))))
871 ;; Split the interval in half
872 (destructuring-bind (y- y
+)
873 (interval-split 0 y t
)
874 (list (interval-expt-< x y-
)
875 (interval-expt-< x y
+))))))
877 ;; The case where x <= 0. Y MUST be an INTEGER for this to work!
878 ;; The calling function must insure this!
879 (loop for interval in
(flatten-list (interval-expt (interval-neg x
) y
))
880 for low
= (interval-low interval
)
881 for high
= (interval-high interval
)
884 collect
(interval-neg interval
)))
886 (destructuring-bind (neg pos
)
887 (interval-split 0 x t t
)
888 (list (interval-expt-< neg y
)
889 (interval-expt-< pos y
))))))
891 ;;; Compute bounds for (expt x y).
892 (defun interval-expt (x y
)
893 (case (interval-range-info x
1)
896 (interval-expt-> x y
))
899 (interval-expt-< x y
))
901 (destructuring-bind (left right
)
902 (interval-split 1 x t t
)
903 (list (interval-expt left y
)
904 (interval-expt right y
))))))
906 (defun fixup-interval-expt (bnd x-int y-int x-type y-type
)
907 (declare (ignore x-int
))
908 ;; Figure out what the return type should be, given the argument
909 ;; types and bounds and the result type and bounds.
910 (cond ((csubtypep x-type
(specifier-type 'integer
))
911 ;; an integer to some power
912 (case (numeric-type-class y-type
)
914 ;; Positive integer to an integer power is either an
915 ;; integer or a rational.
916 (let ((lo (or (interval-low bnd
) '*))
917 (hi (or (interval-high bnd
) '*))
918 (y-lo (interval-low y-int
))
919 (y-hi (interval-high y-int
)))
920 (cond ((and (eq lo
'*)
922 (typep y-lo
'unsigned-byte
)
924 (specifier-type `(integer 0 ,hi
)))
925 ((and (interval-low y-int
)
926 (>= (type-bound-number y-lo
) 0))
928 (specifier-type `(integer ,lo
,hi
)))
930 (specifier-type `(rational ,lo
,hi
))))))
932 ;; Positive integer to rational power is either a rational
933 ;; or a single-float.
934 (let* ((lo (interval-low bnd
))
935 (hi (interval-high bnd
))
937 (floor (type-bound-number lo
))
940 (ceiling (type-bound-number hi
))
942 (f-lo (or (bound-func #'float lo nil
)
944 (f-hi (or (bound-func #'float hi nil
)
946 (specifier-type `(or (rational ,int-lo
,int-hi
)
947 (single-float ,f-lo
, f-hi
)))))
949 ;; A positive integer to a float power is a float.
950 (let ((format (numeric-type-format y-type
)))
952 (modified-numeric-type
954 :low
(coerce-numeric-bound (interval-low bnd
) format
)
955 :high
(coerce-numeric-bound (interval-high bnd
) format
))))
957 ;; A positive integer to a number is a number (for now).
958 (specifier-type 'number
))))
959 ((csubtypep x-type
(specifier-type 'rational
))
960 ;; a rational to some power
961 (case (numeric-type-class y-type
)
963 ;; A positive rational to an integer power is always a rational.
964 (specifier-type `(rational ,(or (interval-low bnd
) '*)
965 ,(or (interval-high bnd
) '*))))
967 ;; A positive rational to rational power is either a rational
968 ;; or a single-float.
969 (let* ((lo (interval-low bnd
))
970 (hi (interval-high bnd
))
972 (floor (type-bound-number lo
))
975 (ceiling (type-bound-number hi
))
977 (f-lo (or (bound-func #'float lo nil
)
979 (f-hi (or (bound-func #'float hi nil
)
981 (specifier-type `(or (rational ,int-lo
,int-hi
)
982 (single-float ,f-lo
, f-hi
)))))
984 ;; A positive rational to a float power is a float.
985 (let ((format (numeric-type-format y-type
)))
987 (modified-numeric-type
989 :low
(coerce-numeric-bound (interval-low bnd
) format
)
990 :high
(coerce-numeric-bound (interval-high bnd
) format
))))
992 ;; A positive rational to a number is a number (for now).
993 (specifier-type 'number
))))
994 ((csubtypep x-type
(specifier-type 'float
))
995 ;; a float to some power
996 (case (numeric-type-class y-type
)
997 ((or integer rational
)
998 ;; A positive float to an integer or rational power is
1000 (let ((format (numeric-type-format x-type
)))
1005 :low
(coerce-numeric-bound (interval-low bnd
) format
)
1006 :high
(coerce-numeric-bound (interval-high bnd
) format
))))
1008 ;; A positive float to a float power is a float of the
1010 (let ((format (float-format-max (numeric-type-format x-type
)
1011 (numeric-type-format y-type
))))
1016 :low
(coerce-numeric-bound (interval-low bnd
) format
)
1017 :high
(coerce-numeric-bound (interval-high bnd
) format
))))
1019 ;; A positive float to a number is a number (for now)
1020 (specifier-type 'number
))))
1022 ;; A number to some power is a number.
1023 (specifier-type 'number
))))
1025 (defun merged-interval-expt (x y
)
1026 (let* ((x-int (numeric-type->interval x
))
1027 (y-int (numeric-type->interval y
)))
1028 (mapcar (lambda (type)
1029 (fixup-interval-expt type x-int y-int x y
))
1030 (flatten-list (interval-expt x-int y-int
)))))
1032 (defun integer-float-p (float)
1034 (multiple-value-bind (significand exponent
) (integer-decode-float float
)
1035 (or (plusp exponent
)
1036 (<= (- exponent
) (sb-kernel::first-bit-set significand
))))))
1038 (defun expt-derive-type-aux (x y same-arg
)
1039 (declare (ignore same-arg
))
1040 (cond ((or (not (numeric-type-real-p x
))
1041 (not (numeric-type-real-p y
)))
1042 ;; Use numeric contagion if either is not real.
1043 (numeric-contagion x y
))
1044 ((or (csubtypep y
(specifier-type 'integer
))
1045 (integer-float-p (nth-value 1 (type-singleton-p y
))))
1046 ;; A real raised to an integer power is well-defined.
1047 (merged-interval-expt x y
))
1048 ;; A real raised to a non-integral power can be a float or a
1050 ((csubtypep x
(specifier-type '(real 0)))
1051 ;; But a positive real to any power is well-defined.
1052 (merged-interval-expt x y
))
1053 ((and (csubtypep x
(specifier-type 'rational
))
1054 (csubtypep y
(specifier-type 'rational
)))
1055 ;; A rational to the power of a rational could be a rational
1056 ;; or a possibly-complex single float
1057 (specifier-type '(or rational single-float
(complex single-float
))))
1059 ;; a real to some power. The result could be a real or a
1061 (float-or-complex-float-type (numeric-contagion x y
)))))
1063 (defoptimizer (expt derive-type
) ((x y
))
1064 (two-arg-derive-type x y
#'expt-derive-type-aux
#'expt
))
1066 ;;; Note we must assume that a type including 0.0 may also include
1067 ;;; -0.0 and thus the result may be complex -infinity + i*pi.
1068 (defun log-derive-type-aux-1 (x)
1069 (elfun-derive-type-simple x
#'log
1070 (if (integer-type-p x
) 0 0d0
)
1072 ;; (log 0) is an error
1073 ;; and there's nothing between 0 and 1 for integers.
1074 (and (integer-type-p x
) 0f0
)
1077 (defun log-derive-type-aux-2 (x y same-arg
)
1078 (let ((log-x (log-derive-type-aux-1 x
))
1079 (log-y (log-derive-type-aux-1 y
))
1080 (accumulated-list nil
))
1081 ;; LOG-X or LOG-Y might be union types. We need to run through
1082 ;; the union types ourselves because /-DERIVE-TYPE-AUX doesn't.
1083 (dolist (x-type (prepare-arg-for-derive-type log-x
))
1084 (dolist (y-type (prepare-arg-for-derive-type log-y
))
1085 (push (/-derive-type-aux x-type y-type same-arg
) accumulated-list
)))
1086 (apply #'type-union
(flatten-list accumulated-list
))))
1088 (defoptimizer (log derive-type
) ((x &optional y
))
1090 (two-arg-derive-type x y
#'log-derive-type-aux-2
#'log
)
1091 (one-arg-derive-type x
#'log-derive-type-aux-1
#'log
)))
1093 (defun atan-derive-type-aux-1 (y)
1094 (elfun-derive-type-simple y
#'atan nil nil
(sb-xc:-
(sb-xc:/ pi
2)) (sb-xc:/ pi
2)))
1096 (defun atan-derive-type-aux-2 (y x same-arg
)
1097 (declare (ignore same-arg
))
1098 ;; The hard case with two args. We just return the max bounds.
1099 (let ((result-type (numeric-contagion y x
)))
1100 (cond ((and (numeric-type-real-p x
)
1101 (numeric-type-real-p y
))
1102 (let* (;; FIXME: This expression for FORMAT seems to
1103 ;; appear multiple times, and should be factored out.
1104 (format (case (numeric-type-class result-type
)
1105 ((integer rational
) 'single-float
)
1106 (t (numeric-type-format result-type
))))
1107 (bound-format (or format
'float
)))
1108 (make-numeric-type :class
'float
1111 :low
(coerce (sb-xc:- pi
) bound-format
)
1112 :high
(coerce pi bound-format
))))
1114 ;; The result is a float or a complex number
1115 (float-or-complex-float-type result-type
)))))
1117 (defoptimizer (atan derive-type
) ((y &optional x
))
1119 (two-arg-derive-type y x
#'atan-derive-type-aux-2
#'atan
)
1120 (one-arg-derive-type y
#'atan-derive-type-aux-1
#'atan
)))
1122 (defun cosh-derive-type-aux (x)
1123 ;; We note that cosh x = cosh |x| for all real x.
1124 (elfun-derive-type-simple
1125 (if (numeric-type-real-p x
)
1126 (abs-derive-type-aux x
)
1128 #'cosh nil nil
0 nil
))
1130 (defoptimizer (cosh derive-type
) ((num))
1131 (one-arg-derive-type num
#'cosh-derive-type-aux
#'cosh
))
1133 (defun phase-derive-type-aux (arg)
1134 (let* ((format (case (numeric-type-class arg
)
1135 ((integer rational
) 'single-float
)
1136 (t (numeric-type-format arg
))))
1137 (bound-type (or format
'float
)))
1138 (cond ((numeric-type-real-p arg
)
1139 (case (interval-range-info> (numeric-type->interval arg
) 0.0)
1141 ;; The number is positive, so the phase is 0.
1142 (make-numeric-type :class
'float
1145 :low
(coerce 0 bound-type
)
1146 :high
(coerce 0 bound-type
)))
1148 ;; The number is always negative, so the phase is pi.
1149 (make-numeric-type :class
'float
1152 :low
(coerce pi bound-type
)
1153 :high
(coerce pi bound-type
)))
1155 ;; We can't tell. The result is 0 or pi. Use a union
1158 (make-numeric-type :class
'float
1161 :low
(coerce 0 bound-type
)
1162 :high
(coerce 0 bound-type
))
1163 (make-numeric-type :class
'float
1166 :low
(coerce pi bound-type
)
1167 :high
(coerce pi bound-type
))))))
1169 ;; We have a complex number. The answer is the range -pi
1170 ;; to pi. (-pi is included because we have -0.)
1171 (make-numeric-type :class
'float
1174 :low
(coerce (sb-xc:- pi
) bound-type
)
1175 :high
(coerce pi bound-type
))))))
1177 (defoptimizer (phase derive-type
) ((num))
1178 (one-arg-derive-type num
#'phase-derive-type-aux
#'phase
))
1180 (deftransform realpart
((x) ((complex rational
)) * :important nil
)
1182 (deftransform imagpart
((x) ((complex rational
)) * :important nil
)
1185 (deftransform realpart
((x) (real) * :important nil
)
1187 (deftransform imagpart
((x) ((and single-float
(not (eql -
0f0
)))) * :important nil
)
1189 (deftransform imagpart
((x) ((and double-float
(not (eql -
0d0
)))) * :important nil
)
1192 ;;; Make REALPART and IMAGPART return the appropriate types. This
1193 ;;; should help a lot in optimized code.
1194 (defun realpart-derive-type-aux (type)
1195 (let ((class (numeric-type-class type
))
1196 (format (numeric-type-format type
)))
1197 (cond ((numeric-type-real-p type
)
1198 ;; The realpart of a real has the same type and range as
1200 (make-numeric-type :class class
1203 :low
(numeric-type-low type
)
1204 :high
(numeric-type-high type
)))
1206 ;; We have a complex number. The result has the same type
1207 ;; as the real part, except that it's real, not complex,
1209 (make-numeric-type :class class
1212 :low
(numeric-type-low type
)
1213 :high
(numeric-type-high type
))))))
1215 (defoptimizer (realpart derive-type
) ((num))
1216 (one-arg-derive-type num
#'realpart-derive-type-aux
#'realpart
))
1218 (defun imagpart-derive-type-aux (type)
1219 (let ((class (numeric-type-class type
))
1220 (format (numeric-type-format type
)))
1221 (cond ((numeric-type-real-p type
)
1222 ;; The imagpart of a real has the same type as the input,
1223 ;; except that it's zero.
1224 (let ((bound-format (or format class
'real
)))
1225 (make-numeric-type :class class
1228 :low
(coerce 0 bound-format
)
1229 :high
(coerce 0 bound-format
))))
1231 ;; We have a complex number. The result has the same type as
1232 ;; the imaginary part, except that it's real, not complex,
1234 (make-numeric-type :class class
1237 :low
(numeric-type-low type
)
1238 :high
(numeric-type-high type
))))))
1240 (defoptimizer (imagpart derive-type
) ((num))
1241 (one-arg-derive-type num
#'imagpart-derive-type-aux
#'imagpart
))
1243 (defun complex-derive-type-aux-1 (re-type)
1244 (if (numeric-type-p re-type
)
1245 (make-numeric-type :class
(numeric-type-class re-type
)
1246 :format
(numeric-type-format re-type
)
1247 :complexp
(if (csubtypep re-type
1248 (specifier-type 'rational
))
1251 :low
(numeric-type-low re-type
)
1252 :high
(numeric-type-high re-type
))
1253 (specifier-type 'complex
)))
1255 (defun complex-derive-type-aux-2 (re-type im-type same-arg
)
1256 (declare (ignore same-arg
))
1257 (if (and (numeric-type-p re-type
)
1258 (numeric-type-p im-type
))
1259 ;; Need to check to make sure numeric-contagion returns the
1260 ;; right type for what we want here.
1262 ;; Also, what about rational canonicalization, like (complex 5 0)
1263 ;; is 5? So, if the result must be complex, we make it so.
1264 ;; If the result might be complex, which happens only if the
1265 ;; arguments are rational, we make it a union type of (or
1266 ;; rational (complex rational)).
1267 (let* ((element-type (numeric-contagion re-type im-type
))
1268 (maybe-rat-result-p (types-equal-or-intersect
1269 element-type
(specifier-type 'rational
)))
1270 (definitely-rat-result-p (csubtypep element-type
(specifier-type 'rational
)))
1271 (real-result-p (and definitely-rat-result-p
1272 (csubtypep im-type
(specifier-type '(eql 0))))))
1274 (real-result-p re-type
)
1276 (let ((complex (specifier-type
1277 `(complex ,(numeric-type-class element-type
)))))
1278 (if (types-equal-or-intersect im-type
(specifier-type '(eql 0)))
1279 (type-union element-type complex
)
1282 (make-numeric-type :class
(numeric-type-class element-type
)
1283 :format
(numeric-type-format element-type
)
1284 :complexp
:complex
))))
1285 (specifier-type 'complex
)))
1287 (defoptimizer (complex derive-type
) ((re &optional im
))
1289 (two-arg-derive-type re im
#'complex-derive-type-aux-2
#'complex
)
1290 (one-arg-derive-type re
#'complex-derive-type-aux-1
#'complex
)))
1292 ;;; Define some transforms for complex operations in lieu of complex operation
1293 ;;; VOPs for most backends. If vops exist, they must support the following
1294 ;;; on complex-single-float and complex-double-float:
1295 ;;; * real-complex, complex-real and complex-complex addition and subtraction
1296 ;;; * complex-real and real-complex multiplication
1297 ;;; * complex-real division
1298 ;;; * sb-vm::swap-complex, which swaps the real and imaginary parts.
1300 ;;; * complex-real, real-complex and complex-complex CL:=
1301 ;;; (complex-complex EQL would usually be a good idea).
1302 (macrolet ((frob (type contagion
)
1304 (deftransform complex
((r) (,type
))
1305 '(complex r
,(coerce 0 type
)))
1306 (deftransform complex
((r i
) (,type
,contagion
))
1307 (when (csubtypep (lvar-type i
) (specifier-type ',type
))
1308 (give-up-ir1-transform))
1309 '(complex r
(truly-the ,type
(coerce i
',type
))))
1310 (deftransform complex
((r i
) (,contagion
,type
))
1311 (when (csubtypep (lvar-type r
) (specifier-type ',type
))
1312 (give-up-ir1-transform))
1313 '(complex (truly-the ,type
(coerce r
',type
)) i
))
1315 ;; Arbitrarily use %NEGATE/COMPLEX-DOUBLE-FLOAT as an indicator
1316 ;; of whether all the operations below are translated by vops.
1317 ;; We could be more fine-grained, but it seems reasonable that
1318 ;; they be implemented on an all-or-none basis.
1319 (unless (vop-existsp :named sb-vm
::%negate
/complex-double-float
)
1321 (deftransform %negate
((z) ((complex ,type
)) * :important nil
)
1322 '(complex (%negate
(realpart z
)) (%negate
(imagpart z
))))
1323 ;; complex addition and subtraction
1324 (deftransform + ((w z
) ((complex ,type
) (complex ,type
)) * :important nil
)
1325 '(complex (+ (realpart w
) (realpart z
))
1326 (+ (imagpart w
) (imagpart z
))))
1327 (deftransform -
((w z
) ((complex ,type
) (complex ,type
)) * :important nil
)
1328 '(complex (- (realpart w
) (realpart z
))
1329 (- (imagpart w
) (imagpart z
))))
1330 ;; Add and subtract a complex and a real.
1331 (deftransform + ((w z
) ((complex ,type
) real
) * :important nil
)
1332 `(complex (+ (realpart w
) z
)
1333 (+ (imagpart w
) ,(coerce 0 ',type
))))
1334 (deftransform + ((z w
) (real (complex ,type
)) * :important nil
)
1335 `(complex (+ (realpart w
) z
)
1336 (+ (imagpart w
) ,(coerce 0 ',type
))))
1337 ;; Add and subtract a real and a complex number.
1338 (deftransform -
((w z
) ((complex ,type
) real
) * :important nil
)
1339 `(complex (- (realpart w
) z
)
1340 (- (imagpart w
) ,(coerce 0 ',type
))))
1341 (deftransform -
((z w
) (real (complex ,type
)) * :important nil
)
1342 `(complex (- z
(realpart w
))
1343 (- ,(coerce 0 ',type
) (imagpart w
))))
1344 ;; Multiply a complex by a real or vice versa.
1345 (deftransform * ((w z
) ((complex ,type
) real
) * :important nil
)
1346 '(complex (* (realpart w
) z
) (* (imagpart w
) z
)))
1347 (deftransform * ((z w
) (real (complex ,type
)) * :important nil
)
1348 '(complex (* (realpart w
) z
) (* (imagpart w
) z
)))
1349 ;; conjugate of complex number
1350 (deftransform conjugate
((z) ((complex ,type
)) * :important nil
)
1351 '(complex (realpart z
) (- (imagpart z
))))
1353 (deftransform = ((w z
) ((complex ,type
) (complex ,type
)) * :important nil
)
1354 '(and (= (realpart w
) (realpart z
))
1355 (= (imagpart w
) (imagpart z
))))
1356 (deftransform = ((w z
) ((complex ,type
) real
) * :important nil
)
1357 '(and (= (realpart w
) z
) (zerop (imagpart w
))))
1358 (deftransform = ((w z
) (real (complex ,type
)) * :important nil
)
1359 '(and (= (realpart z
) w
) (zerop (imagpart z
))))
1360 ;; Multiply two complex numbers.
1361 (deftransform * ((x y
) ((complex ,type
) (complex ,type
)) * :important nil
)
1362 '(let* ((rx (realpart x
))
1366 (complex (- (* rx ry
) (* ix iy
))
1367 (+ (* rx iy
) (* ix ry
)))))
1368 ;; Divide a complex by a real.
1369 (deftransform / ((w z
) ((complex ,type
) real
) * :important nil
)
1370 '(complex (/ (realpart w
) z
) (/ (imagpart w
) z
)))
1373 ;; Divide two complex numbers.
1374 (deftransform / ((x y
) ((complex ,type
) (complex ,type
)) * :important nil
)
1375 (if (vop-existsp :translate sb-vm
::swap-complex
)
1376 '(let* ((cs (conjugate (sb-vm::swap-complex x
)))
1379 (if (> (abs ry
) (abs iy
))
1380 (let* ((r (/ iy ry
))
1381 (dn (+ ry
(* r iy
))))
1382 (/ (+ x
(* cs r
)) dn
))
1383 (let* ((r (/ ry iy
))
1384 (dn (+ iy
(* r ry
))))
1385 (/ (+ (* x r
) cs
) dn
))))
1386 '(let* ((rx (realpart x
))
1390 (if (> (abs ry
) (abs iy
))
1391 (let* ((r (/ iy ry
))
1392 (dn (+ ry
(* r iy
))))
1393 (complex (/ (+ rx
(* ix r
)) dn
)
1394 (/ (- ix
(* rx r
)) dn
)))
1395 (let* ((r (/ ry iy
))
1396 (dn (+ iy
(* r ry
))))
1397 (complex (/ (+ (* rx r
) ix
) dn
)
1398 (/ (- (* ix r
) rx
) dn
)))))))
1399 ;; Divide a real by a complex.
1400 (deftransform / ((x y
) (real (complex ,type
)) * :important nil
)
1401 (if (vop-existsp :translate sb-vm
::swap-complex
)
1402 '(let* ((ry (realpart y
))
1404 (if (> (abs ry
) (abs iy
))
1405 (let* ((r (/ iy ry
))
1406 (dn (+ ry
(* r iy
))))
1407 (/ (complex x
(- (* x r
))) dn
))
1408 (let* ((r (/ ry iy
))
1409 (dn (+ iy
(* r ry
))))
1410 (/ (complex (* x r
) (- x
)) dn
))))
1411 '(let* ((ry (realpart y
))
1413 (if (> (abs ry
) (abs iy
))
1414 (let* ((r (/ iy ry
))
1415 (dn (+ ry
(* r iy
))))
1417 (/ (- (* x r
)) dn
)))
1418 (let* ((r (/ ry iy
))
1419 (dn (+ iy
(* r ry
))))
1420 (complex (/ (* x r
) dn
)
1423 (deftransform cis
((z) ((,type
)) *)
1424 '(complex (cos z
) (sin z
)))
1426 (frob single-float
(or rational single-float
))
1427 (frob double-float
(or rational single-float double-float
)))
1430 ;;;; float contagion
1431 (deftransform single-float-real-contagion
((x y
) * * :node node
:defun-only t
)
1432 (if (csubtypep (lvar-type y
) (specifier-type 'single-float
))
1433 (give-up-ir1-transform)
1434 `(,(lvar-fun-name (basic-combination-fun node
)) x
(%single-float y
))))
1436 (deftransform real-single-float-contagion
((x y
) * * :node node
:defun-only t
)
1437 (if (csubtypep (lvar-type x
) (specifier-type 'single-float
))
1438 (give-up-ir1-transform)
1439 `(,(lvar-fun-name (basic-combination-fun node
)) (%single-float x
) y
)))
1441 (deftransform double-float-real-contagion
((x y
) * * :node node
:defun-only t
)
1442 (if (csubtypep (lvar-type y
) (specifier-type 'double-float
))
1443 (give-up-ir1-transform)
1444 `(,(lvar-fun-name (basic-combination-fun node
)) x
(%double-float y
))))
1446 (deftransform real-double-float-contagion
((x y
) * * :node node
:defun-only t
)
1447 (if (csubtypep (lvar-type x
) (specifier-type 'double-float
))
1448 (give-up-ir1-transform)
1449 `(,(lvar-fun-name (basic-combination-fun node
)) (%double-float x
) y
)))
1451 (deftransform double-float-real-contagion-cmp
((x y
) * * :node node
:defun-only t
)
1452 (cond ((csubtypep (lvar-type y
) (specifier-type 'double-float
))
1453 (give-up-ir1-transform))
1454 ;; Turn (= single-float 1d0) into (= single-float 1f0)
1455 ((and (constant-lvar-p x
)
1456 (csubtypep (lvar-type y
) (specifier-type 'single-float
))
1457 (let ((x (lvar-value x
)))
1458 (when (and (safe-single-coercion-p x
)
1459 (= x
(coerce x
'single-float
)))
1460 `(,(lvar-fun-name (basic-combination-fun node
)) ,(coerce x
'single-float
) y
)))))
1462 `(,(lvar-fun-name (basic-combination-fun node
)) x
(%double-float y
)))))
1464 (deftransform real-double-float-contagion-cmp
((x y
) * * :node node
:defun-only t
)
1465 (cond ((csubtypep (lvar-type x
) (specifier-type 'double-float
))
1466 (give-up-ir1-transform))
1467 ((and (constant-lvar-p y
)
1468 (csubtypep (lvar-type x
) (specifier-type 'single-float
))
1469 (let ((y (lvar-value y
)))
1470 (when (and (safe-single-coercion-p y
)
1471 (= y
(coerce y
'single-float
)))
1472 `(,(lvar-fun-name (basic-combination-fun node
)) x
,(coerce y
'single-float
))))))
1474 `(,(lvar-fun-name (basic-combination-fun node
)) (%double-float x
) y
))))
1477 (%deftransform op nil
'(function (single-float real
) single-float
)
1478 #'single-float-real-contagion nil
)
1479 (%deftransform op nil
'(function (real single-float
) single-float
)
1480 #'real-single-float-contagion nil
)
1481 (%deftransform op nil
'(function (double-float real
))
1482 #'double-float-real-contagion nil
)
1483 (%deftransform op nil
'(function (real double-float
))
1484 #'real-double-float-contagion nil
)
1486 (%deftransform op nil
'(function ((complex single-float
) real
) (complex single-float
))
1487 #'single-float-real-contagion nil
)
1488 (%deftransform op nil
'(function (real (complex single-float
)) (complex single-float
))
1489 #'real-single-float-contagion nil
)
1490 (%deftransform op nil
'(function ((complex double-float
) real
) (complex double-float
))
1491 #'double-float-real-contagion nil
)
1492 (%deftransform op nil
'(function (real (complex double-float
)) (complex double-float
))
1493 #'real-double-float-contagion nil
)))
1494 (dolist (op '(+ * / -
))
1498 (%deftransform op nil
'(function (single-float real
) (values t single-float
))
1499 #'single-float-real-contagion nil
)
1500 (%deftransform op nil
'(function (real single-float
) (values t single-float
))
1501 #'real-single-float-contagion nil
)
1502 (%deftransform op nil
'(function (double-float real
))
1503 #'double-float-real-contagion nil
)
1504 (%deftransform op nil
'(function (real double-float
))
1505 #'real-double-float-contagion nil
)))
1506 (dolist (op '(floor ceiling round truncate ffloor fceiling fround ftruncate
))
1510 (%deftransform op nil
`(function (single-float (integer ,most-negative-exactly-single-float-integer
1511 ,most-positive-exactly-single-float-integer
)))
1512 #'single-float-real-contagion nil
)
1513 (%deftransform op nil
`(function ((integer ,most-negative-exactly-single-float-integer
1514 ,most-positive-exactly-single-float-integer
)
1516 #'real-single-float-contagion nil
)
1518 (%deftransform op nil
`(function (double-float
1520 (integer ,most-negative-exactly-double-float-integer
1521 ,most-positive-exactly-double-float-integer
))))
1522 #'double-float-real-contagion-cmp nil
)
1523 (%deftransform op nil
`(function ((or single-float
1524 (integer ,most-negative-exactly-double-float-integer
1525 ,most-positive-exactly-double-float-integer
))
1527 #'real-double-float-contagion-cmp nil
)))
1528 (dolist (op '(= < > <= >=))
1531 (%deftransform
'= nil
'(function ((complex double-float
) single-float
))
1532 #'double-float-real-contagion nil
)
1533 (%deftransform
'= nil
'(function (single-float (complex double-float
)))
1534 #'real-double-float-contagion nil
)
1536 (deftransform complex
((realpart &optional imagpart
) (rational &optional
(or null
(integer 0 0))) * :important nil
)
1539 ;;; Here are simple optimizers for SIN, COS, and TAN. They do not
1540 ;;; produce a minimal range for the result; the result is the widest
1541 ;;; possible answer. This gets around the problem of doing range
1542 ;;; reduction correctly but still provides useful results when the
1543 ;;; inputs are union types.
1544 (defun trig-derive-type-aux (arg domain fun
1545 &optional def-lo def-hi
(increasingp t
))
1548 (flet ((floatify-format ()
1549 (case (numeric-type-class arg
)
1550 ((integer rational
) 'single-float
)
1551 (t (numeric-type-format arg
)))))
1552 (cond ((eq (numeric-type-complexp arg
) :complex
)
1553 (make-numeric-type :class
'float
1554 :format
(floatify-format)
1555 :complexp
:complex
))
1556 ((numeric-type-real-p arg
)
1557 (let* ((format (floatify-format))
1558 (bound-type (or format
'float
)))
1559 ;; If the argument is a subset of the "principal" domain
1560 ;; of the function, we can compute the bounds because
1561 ;; the function is monotonic. We can't do this in
1562 ;; general for these periodic functions because we can't
1563 ;; (and don't want to) do the argument reduction in
1564 ;; exactly the same way as the functions themselves do
1566 (if (csubtypep arg domain
)
1567 (let ((res-lo (bound-func fun
(numeric-type-low arg
) nil
))
1568 (res-hi (bound-func fun
(numeric-type-high arg
) nil
)))
1570 (rotatef res-lo res-hi
))
1574 :low
(coerce-numeric-bound res-lo bound-type
)
1575 :high
(coerce-numeric-bound res-hi bound-type
)))
1579 :low
(and def-lo
(coerce def-lo bound-type
))
1580 :high
(and def-hi
(coerce def-hi bound-type
))))))
1582 (float-or-complex-float-type arg def-lo def-hi
)))))))
1584 (defoptimizer (sin derive-type
) ((num))
1585 (one-arg-derive-type
1588 ;; Derive the bounds if the arg is in [-pi/2, pi/2].
1589 (trig-derive-type-aux
1591 (specifier-type `(float ,(sb-xc:-
(sb-xc:/ pi
2)) ,(sb-xc:/ pi
2)))
1596 (defoptimizer (cos derive-type
) ((num))
1597 (one-arg-derive-type
1600 ;; Derive the bounds if the arg is in [0, pi].
1601 (trig-derive-type-aux arg
1602 (specifier-type `(float 0d0
,pi
))
1608 (defoptimizer (tan derive-type
) ((num))
1609 (one-arg-derive-type
1612 ;; Derive the bounds if the arg is in [-pi/2, pi/2].
1613 (trig-derive-type-aux arg
1614 (specifier-type `(float ,(sb-xc:-
(sb-xc:/ pi
2))
1620 (defoptimizer (conjugate derive-type
) ((num))
1621 (one-arg-derive-type num
1623 (flet ((most-negative-bound (l h
)
1625 (if (< (type-bound-number l
) (- (type-bound-number h
)))
1627 (set-bound (- (type-bound-number h
)) (consp h
)))))
1628 (most-positive-bound (l h
)
1630 (if (> (type-bound-number h
) (- (type-bound-number l
)))
1632 (set-bound (- (type-bound-number l
)) (consp l
))))))
1633 (if (numeric-type-real-p arg
)
1635 (let ((low (numeric-type-low arg
))
1636 (high (numeric-type-high arg
)))
1637 (let ((new-low (most-negative-bound low high
))
1638 (new-high (most-positive-bound low high
)))
1639 (modified-numeric-type arg
:low new-low
:high new-high
))))))
1642 (defoptimizer (cis derive-type
) ((num))
1643 (one-arg-derive-type num
1646 `(complex ,(or (numeric-type-format arg
) 'float
))))
1650 ;;;; TRUNCATE, FLOOR, CEILING, and ROUND
1651 (deftransform truncate
((x &optional by
)
1652 (t &optional
(constant-arg (member 1))))
1653 '(unary-truncate x
))
1655 (deftransform round
((x &optional by
)
1656 (t &optional
(constant-arg (member 1))))
1657 '(let ((res (%unary-round x
)))
1658 (values res
(locally
1659 (declare (flushable %single-float
1663 (deftransform %unary-truncate
((x) (single-float))
1664 `(values (unary-truncate x
)))
1665 (deftransform %unary-truncate
((x) (double-float))
1666 `(values (unary-truncate x
)))
1668 (defun value-within-numeric-type (type)
1670 (when (ctypep x type
)
1671 (return-from value-within-numeric-type x
)))
1673 (multiple-value-bind (frac exp sign
)
1674 (integer-decode-float float
)
1675 (* (scale-float (float (1+ frac
) float
) exp
)
1678 (multiple-value-bind (frac exp sign
)
1679 (integer-decode-float float
)
1680 (* (scale-float (float (1- frac
) float
) exp
)
1698 (ratio-between (low high
)
1699 (+ low
(/ (- high low
) 2)))
1701 (when (numeric-type-p x
)
1702 (let ((lo (numeric-type-low x
))
1703 (hi (numeric-type-high x
)))
1709 (try (next (car lo
))))
1711 (try (prev (car hi
))))
1712 (when (and (typep lo
'(cons rational
))
1713 (typep hi
'(cons rational
)))
1714 (try (ratio-between (car lo
) (car hi
))))
1715 (when (csubtypep x
(specifier-type 'rational
))
1717 (when (csubtypep x
(specifier-type 'double-float
))
1719 (when (csubtypep x
(specifier-type 'single-float
))
1722 (numeric-type (numeric type
))
1723 (union-type (mapc #'numeric
(union-type-types type
))))
1724 (error "Couldn't come up with a value for ~s" type
)))
1726 #-
(or sb-xc-host
64-bit
)
1728 (declaim (inline %make-double-float
))
1729 (defun %make-double-float
(bits)
1730 (make-double-float (ash bits -
32) (ldb (byte 32 0) bits
))))
1732 ;;; Transform inclusive integer bounds so that they work on floats
1733 ;;; before truncating to zero.
1736 `(defun ,(symbolicate type
'-integer-bounds
) (low high
)
1737 (macrolet ((const (name)
1738 (package-symbolicate :sb-vm
',type
'- name
)))
1739 (labels ((fractions-p (number)
1740 (< (integer-length (abs number
))
1744 (sb-kernel:make-single-float
0)
1745 (let* ((negative (minusp number
))
1746 (number (abs number
))
1747 (length (integer-length number
))
1748 (shift (- length
(const digits
)))
1749 (shifted (truly-the fixnum
1752 ;; Cut off the hidden bit
1753 (signif (ldb (const significand-byte
) shifted
))
1754 (exp (+ (const bias
) length
))
1756 (byte-position (const exponent-byte
)))))
1758 ;; If rounding up overflows this will increase the exponent too
1759 (let ((bits (+ bits signif
)))
1761 (setf bits
(logior (ash -
1 ,(case type
1766 (single-float 'make-single-float
)
1767 (double-float '%make-double-float
)) bits
))))))
1768 (values (if (<= low
0)
1769 (if (fractions-p low
)
1775 (if (fractions-p high
)
1781 (deftransform unary-truncate
((x) * * :result result
:node node
)
1782 (delay-ir1-transform node
:constraint
)
1783 (unless (or (lvar-single-value-p result
)
1784 (mv-bind-unused-p result
1))
1785 (give-up-ir1-transform))
1786 (let ((rem-type (second (values-type-required (node-derived-type node
)))))
1787 `(values (%unary-truncate x
)
1788 ,(value-within-numeric-type rem-type
))))
1790 (macrolet ((def (type)
1791 `(deftransform unary-truncate
((number) (,type
) * :node node
)
1792 (let ((cast (cast-or-check-bound-type node
)))
1794 (csubtypep cast
(specifier-type 'sb-vm
:signed-word
)))
1795 (let ((int (type-approximate-interval cast
)))
1797 (multiple-value-bind (low high
) (,(symbolicate type
'-integer-bounds
)
1799 (interval-high int
))
1801 '(,',type
,low
,high
))
1802 (let ((truncated (truly-the ,(type-specifier cast
) (,',(symbolicate '%unary-truncate
/ type
) number
))))
1803 (declare (flushable ,',(symbolicate "%" type
)))
1806 (coerce truncated
',',type
))))
1807 ,(internal-type-error-call 'number
(type-specifier cast
) 'truncate-to-integer
)))))
1810 ,(symbol-value (package-symbolicate :sb-kernel
'most-negative-fixnum- type
))
1811 ,(symbol-value (package-symbolicate :sb-kernel
'most-positive-fixnum- type
))))
1812 (let ((truncated (truly-the fixnum
(,(symbolicate '%unary-truncate
/ type
) number
))))
1813 (declare (flushable ,(symbolicate "%" type
)))
1816 (coerce truncated
',type
))))
1817 (,(symbolicate 'unary-truncate- type
'-to-bignum
) number
)))))))
1822 (macrolet ((def (type other-float-arg-types
)
1823 (let* ((unary (symbolicate "%UNARY-TRUNCATE/" type
))
1824 (unary-to-bignum (symbolicate '%unary-truncate- type
'-to-bignum
))
1825 (coerce (symbolicate "%" type
))
1826 (unary `(lambda (number)
1829 ,(symbol-value (package-symbolicate :sb-kernel
'most-negative-fixnum- type
))
1830 ,(symbol-value (package-symbolicate :sb-kernel
'most-positive-fixnum- type
))))
1831 (let ((r (truly-the fixnum
(,unary number
))))
1833 (declare (flushable ,coerce
))
1835 (let ((r (,unary-to-bignum number
)))
1837 number
;; no fractional part
1840 (declare (flushable ,coerce
))
1842 `(deftransform truncate
((x &optional y
)
1844 &optional
(or ,type
,@other-float-arg-types integer
))
1845 * :result result
:node node
)
1846 (let* ((result-type (and result
1847 (lvar-derived-type result
)))
1848 (compute-all (and (or (eq result-type
*wild-type
*)
1849 (values-type-p result-type
))
1850 (not (type-single-value-p result-type
))))
1852 (and (constant-lvar-p y
) (sb-xc:= 1 (lvar-value y
))))))
1856 `(let ((res (,',unary x
)))
1857 ;; Dummy secondary value!
1860 `(let* ((f (,',coerce y
))
1862 (multiple-value-bind (res float-res
) (,',unary div
)
1864 (- x
(* f float-res
)))))
1865 `(let* ((f (,',coerce y
))
1866 (res (,',unary
(/ x f
))))
1867 ;; Dummy secondary value!
1868 (values res x
)))))))))
1869 (def single-float
())
1870 (def double-float
(single-float)))
1872 (macrolet ((def (name type other-float-arg-types
)
1873 (let* ((unary-to-bignum (symbolicate 'unary-truncate- type
'-to-bignum
))
1874 (coerce (symbolicate "%" type
))
1875 (fixnum-type `(,type
1876 ,(symbol-value (package-symbolicate :sb-kernel
'most-negative-fixnum- type
))
1877 ,(symbol-value (package-symbolicate :sb-kernel
'most-positive-fixnum- type
)))))
1878 `(deftransform ,name
((number &optional divisor
)
1880 &optional
(or ,type
,@other-float-arg-types integer
)))
1881 (let ((one-p (or (not divisor
)
1882 (and (constant-lvar-p divisor
) (sb-xc:= (lvar-value divisor
) 1)))))
1885 `((f-divisor (,',coerce divisor
))
1886 (div (/ number f-divisor
))))
1889 (double-float 'round-double
)
1890 (single-float 'round-single
))
1891 div
,,(keywordicate name
))))
1892 (values (if (typep div
',',fixnum-type
)
1893 ,',(if-vop-existsp (:translate %unary-ceiling
)
1894 `(truly-the fixnum
(,(symbolicate '%unary- name
) div
))
1895 `(%unary-truncate
(truly-the ,fixnum-type quot
)))
1896 (,',unary-to-bignum quot
))
1897 (- number
(* ,@(unless one-p
1902 (double-float 0.0d0
)
1903 (single-float 0.0f0
))))))))))))
1904 (def floor single-float
())
1905 (def floor double-float
(single-float))
1906 (def ceiling single-float
())
1907 (def ceiling double-float
(single-float))
1909 (def truncate single-float
())
1911 (def truncate double-float
(single-float)))
1915 (defknown (%unary-fround
) (real) float
(movable foldable flushable
))
1917 (defknown (%unary-fround
/double
) (double-float) double-float
1918 (movable foldable flushable
))
1920 (deftransform %unary-fround
((x) (single-float))
1921 `(cond ((typep x
'(or (single-float -
0.5f0
(0f0)) (eql -
0f0
)))
1923 ((typep x
'(single-float ,(float (- (expt 2 sb-vm
:single-float-digits
)) 1f0
)
1924 ,(float (1- (expt 2 sb-vm
:single-float-digits
)) 1f0
)))
1925 (float (round x
) 1f0
))
1930 (deftransform %unary-fround
((x) (double-float))
1931 `(cond ((typep x
'(or (double-float -
0.5d0
(0d0)) (eql -
0d0
)))
1933 ((typep x
'(double-float ,(float (- (expt 2 sb-vm
:double-float-digits
)) 1d0
)
1934 ,(float (1- (expt 2 sb-vm
:double-float-digits
)) 1d0
)))
1935 (float (round x
) 1d0
))
1942 (defun %unary-fround
/double
(x)
1943 (declare (muffle-conditions compiler-note
))
1944 (declare (type double-float x
))
1945 (declare (optimize speed
(safety 0)))
1946 (let* ((high (double-float-high-bits x
))
1947 (low (double-float-low-bits x
))
1948 (exp (ldb sb-vm
:double-float-hi-exponent-byte high
))
1949 (biased (the double-float-exponent
1950 (- exp sb-vm
:double-float-bias
))))
1951 (declare (type (signed-byte 32) high
)
1952 (type (unsigned-byte 32) low
))
1954 ((= exp sb-vm
:double-float-normal-exponent-max
) x
)
1955 ((<= biased -
1) (* x
0d0
)) ; [0,0.5)
1956 ((and (= biased
0) (= low
0) (= (ldb sb-vm
:double-float-hi-significand-byte high
) 0)) ; [0.5,0.5]
1958 ((= biased
0) (float-sign x
1d0
)) ; (0.5,1.0)
1959 ((= biased
1) ; [1.0,2.0)
1961 ((>= (ldb sb-vm
:double-float-hi-significand-byte high
) (ash 1 19))
1963 (t (float-sign x
1d0
))))
1964 ((>= biased
(float-digits x
)) x
)
1966 ;; it's probably possible to do something very contorted
1967 ;; to avoid consing intermediate bignums, by performing
1968 ;; arithmetic on the fractional part, the low integer
1969 ;; part, the high integer part, and the exponent of the
1970 ;; double float. But in the interest of getting
1971 ;; something correct to start with, delegate to ROUND.
1972 (float (round x
) 1d0
)))))
1973 (deftransform %unary-fround
((x) (double-float))
1974 `(%unary-fround
/double x
)))
1976 (macrolet ((def (name type
&optional
(suffix ""))
1977 `(deftransform ,name
((x mode
) (t (constant-arg t
)))
1978 (let ((fun (case (lvar-value mode
)
1979 (:floor
,(format nil
"floor~a" suffix
))
1980 (:ceiling
,(format nil
"ceil~a" suffix
))
1981 (:truncate
,(format nil
"trunc~a" suffix
)))))
1983 (declare (optimize (sb-c:alien-funcall-saves-fp-and-pc
0)))
1984 (alien-funcall (%alien-value
1985 (foreign-symbol-sap ,fun nil
) 0
1986 ,(parse-alien-type '(function ,type
,type
) nil
))
1988 (def round-single single-float
"f")
1989 (def round-double double-float
)))
1992 (deftransform fround
((number &optional divisor
) (double-float &optional t
))
1993 (if (or (not divisor
)
1994 (and (constant-lvar-p divisor
)
1995 (= (lvar-value divisor
) 1)))
1996 `(let ((res (round-double number
:round
)))
1997 (values res
(- number res
)))
1998 `(let* ((divisor (%double-float divisor
))
1999 (res (round-double (/ number
(%double-float divisor
)) :round
)))
2000 (values res
(- number
(* res divisor
))))))
2003 (deftransform fround
((number &optional divisor
) (single-float &optional
(or null single-float rational
)))
2004 (if (or (not divisor
)
2005 (and (constant-lvar-p divisor
)
2006 (= (lvar-value divisor
) 1)))
2007 `(let ((res (round-single number
:round
)))
2008 (values res
(- number res
)))
2009 `(let* ((divisor (%single-float divisor
))
2010 (res (round-single (/ number divisor
) :round
)))
2011 (values res
(- number
(* res divisor
))))))
2015 ;;; Dumping of double-float literals in genesis got some bits messed up,
2016 ;;; but only if the double-float was the value of a slot in a ctype instance.
2017 ;;; It was broken for either endianness, but miraculously didn't crash
2018 ;;; for little-endian builds even though it could have.
2019 ;;; (The dumped constants were legal normalalized float bit patterns, albeit wrong)
2020 ;;; For 32-bit big-endian machines, the bit patterns were those of subnormals.
2021 ;;; So thank goodness for that - it allowed detection of the problem.
2022 (defun test-ctype-involving-double-float ()
2023 (specifier-type '(double-float #.pi
)))
2024 (assert (sb-xc:= (numeric-type-low (test-ctype-involving-double-float)) pi
))
2026 ;;; Dummy functions to test that complex number are dumped correctly in genesis.
2027 (defun try-folding-complex-single ()
2028 (let ((re (make-single-float #x4E000000
))
2029 (im (make-single-float #x-21800000
)))
2030 (values (complex re im
)
2031 (locally (declare (notinline complex
)) (complex re im
)))))
2033 (defun try-folding-complex-double ()
2034 (let ((re (make-double-float #X3FE62E42
#xFEFA39EF
))
2035 (im (make-double-float #X43CFFFFF
#XFFFFFFFF
)))
2036 (values (complex re im
)
2037 (locally (declare (notinline complex
)) (complex re im
)))))
2039 (dolist (test '(try-folding-complex-single try-folding-complex-double
))
2040 (multiple-value-bind (a b
) (funcall test
)
2043 (let ((code (fun-code-header (symbol-function test
))))
2044 (aver (loop for index from sb-vm
:code-constants-offset
2045 below
(code-header-words code
)
2046 thereis
(typep (code-header-ref code index
) 'complex
))))
2049 (defun more-folding ()
2050 (values (complex single-float-positive-infinity single-float-positive-infinity
)
2051 (complex single-float-negative-infinity single-float-positive-infinity
)
2052 (complex single-float-negative-infinity single-float-negative-infinity
)
2053 (complex single-float-positive-infinity single-float-negative-infinity
)))
2055 (multiple-value-bind (a b c d
) (funcall 'more-folding
)
2056 (assert (sb-ext:float-infinity-p
(realpart a
)))
2057 (assert (sb-ext:float-infinity-p
(imagpart a
)))
2058 (assert (sb-ext:float-infinity-p
(realpart b
)))
2059 (assert (sb-ext:float-infinity-p
(imagpart b
)))
2060 (assert (sb-ext:float-infinity-p
(realpart c
)))
2061 (assert (sb-ext:float-infinity-p
(imagpart c
)))
2062 (assert (sb-ext:float-infinity-p
(realpart d
)))
2063 (assert (sb-ext:float-infinity-p
(imagpart d
)))
2065 (let ((code (fun-code-header (symbol-function 'more-folding
))))
2066 (aver (loop for index from sb-vm
:code-constants-offset
2067 below
(code-header-words code
)
2068 thereis
(typep (code-header-ref code index
) 'complex
))))
2069 (fmakunbound 'more-folding
))
2071 ;;; Inline (= float 1) by doing two comparisons.
2072 (macrolet ((def (op)
2073 `(deftransform ,op
((x y
) (:or
((float
2074 (integer #.most-negative-exactly-single-float-integer
2075 #.most-positive-exactly-single-float-integer
)) *)
2076 ((float (constant-arg float
)) *))
2077 * :node node
:important nil
2078 :policy
(> speed
1))
2079 (unless (and (types-equal-or-intersect (lvar-type x
) (specifier-type 'double-float
))
2080 (types-equal-or-intersect (lvar-type x
) (specifier-type 'single-float
)))
2081 (give-up-ir1-transform))
2082 (delay-ir1-transform node
:ir1-phases
)
2083 (if (csubtypep (lvar-type y
) (specifier-type 'float
))
2084 (let ((y (lvar-value y
)))
2085 (if (and (safe-single-coercion-p y
)
2086 (sb-xc:= y
(coerce y
'single-float
))
2087 (sb-xc:= y
(coerce y
'double-float
)))
2088 `(if (single-float-p x
)
2089 (,',op
(truly-the single-float x
) ,(coerce y
'single-float
))
2090 (,',op
(truly-the double-float x
) ,(coerce y
'double-float
)))
2091 (give-up-ir1-transform)))
2092 `(if (single-float-p x
)
2093 (,',op
(truly-the single-float x
) (%single-float y
))
2094 (,',op
(truly-the double-float x
) (%double-float y
)))))))
2101 (deftransform phase
((n))
2102 (splice-fun-args n
'complex
2)
2106 (defoptimizer (atan externally-checkable-type
) ((y &rest x
) node
)
2108 (specifier-type 'real
)
2109 (specifier-type 'number
)))
