Newer
Older
;;; -*- Mode: Lisp; Package: C; Log: code.log -*-
;;;
;;; **********************************************************************
;;; This code was written as part of the CMU Common Lisp project at
;;; Carnegie Mellon University, and has been placed in the public domain.
;;;
(ext:file-comment
"$Header: src/compiler/float-tran.lisp $")
;;; **********************************************************************
;;;
;;; This file contains floating-point specific transforms, and may be somewhat
;;; implementation dependent in its assumptions of what the formats are.
;;;
;;; Author: Rob MacLachlan
;;;
(in-package "C")
(intl:textdomain "cmucl")
;;;; Coercions:
(defknown %single-float (real) single-float (movable foldable flushable))
(defknown %double-float (real) double-float (movable foldable flushable))
(deftransform float ((n prototype) (* single-float) * :when :both)
'(%single-float n))
(deftransform float ((n prototype) (* double-float) * :when :both)
(deftransform float ((n) *)
`(if (floatp n) n (%single-float n)))
(deftransform %single-float ((n) (single-float) * :when :both)
(deftransform %double-float ((n) (double-float) * :when :both)
#+double-double
(progn
(defknown %double-double-float (real)
double-double-float
(movable foldable flushable))
(deftransform float ((n prototype) (* double-double-float) * :when :both)
'(%double-double-float n))
(deftransform %double-float ((n) (double-double-float) * :when :both)
'(double-double-hi n))
(deftransform %single-float ((n) (double-double-float) * :when :both)
'(float (double-double-hi n) 1f0))
(deftransform %double-double-float ((n) (double-double-float) * :when :both)
'n)
#+nil
(defun %double-double-float (n)
(make-double-double-float (float n 1d0) 0d0))
;; Moved to code/float.lisp, because we need this relatively early in
;; the build process to handle float and real types.
#+nil
(defun %double-double-float (n)
(typecase n
(fixnum
(%make-double-double-float (float n 1d0) 0d0))
(%make-double-double-float (float n 1d0) 0d0))
(%make-double-double-float (float n 1d0) 0d0))
(double-double-float
n)
(bignum
(bignum:bignum-to-float n 'double-double-float))
(ratio
(kernel::float-ratio n 'double-double-float))))
); progn
(defknown %complex-single-float (number) (complex single-float)
(movable foldable flushable))
(defknown %complex-double-float (number) (complex double-float)
(movable foldable flushable))
(defknown %complex-double-double-float (number) (complex double-double-float)
(movable foldable flushable))
(macrolet
((frob (type)
(let ((name (symbolicate "%COMPLEX-" type "-FLOAT"))
(convert (symbolicate "%" type "-FLOAT")))
`(progn
(defun ,name (n)
(declare (number n))
(etypecase n
(real
(complex (,convert n)))
(complex
(complex (,convert (realpart n))
(,convert (imagpart n))))))
(deftransform ,name ((n) ((complex ,(symbolicate type "-FLOAT"))) * :when :both)
'n)))))
(frob single)
(frob double)
#+double-double
(frob double-double))
(deftransform coerce ((n type) (* *) * :when :both)
(unless (constant-continuation-p type)
(give-up))
`(the ,(continuation-value type)
,(let ( (tspec (specifier-type (continuation-value type))) )
(cond #+double-double
((csubtypep tspec (specifier-type 'double-double-float))
'(%double-double-float n))
((csubtypep tspec (specifier-type 'double-float))
'(%double-float n))
((csubtypep tspec (specifier-type 'float))
'(%single-float n))
#+double-double
((csubtypep tspec (specifier-type '(complex double-double-float)))
'(%complex-double-double-float n))
((csubtypep tspec (specifier-type '(complex double-float)))
'(%complex-double-float n))
((csubtypep tspec (specifier-type '(complex single-float)))
'(%complex-single-float n))
;;; Not strictly float functions, but primarily useful on floats:
(macrolet ((frob (fun ufun)
`(progn
(defknown ,ufun (real) integer (movable foldable flushable))
(deftransform ,fun ((x &optional by)
(* &optional
(constant-argument (member 1))))
'(let ((res (,ufun x)))
(values res (- x res)))))))
(frob round %unary-round))
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
(defknown %unary-truncate (real) integer
(movable foldable flushable))
;; Convert (truncate x y) to the obvious implementation. We only want
;; this when under certain conditions and let the generic truncate
;; handle the rest. (Note: if y = 1, the divide and multiply by y
;; should be removed by other deftransforms.)
(deftransform truncate ((x &optional y)
(float &optional (or float integer)))
'(let ((res (%unary-truncate (/ x y))))
(values res (- x (* y res)))))
(deftransform floor ((number &optional divisor)
(float &optional (or integer float)))
'(multiple-value-bind (tru rem) (truncate number divisor)
(if (and (not (zerop rem))
(if (minusp divisor)
(plusp number)
(minusp number)))
(values (1- tru) (+ rem divisor))
(values tru rem))))
(deftransform ceiling ((number &optional divisor)
(float &optional (or integer float)))
'(multiple-value-bind (tru rem) (truncate number divisor)
(if (and (not (zerop rem))
(if (minusp divisor)
(minusp number)
(plusp number)))
(values (1+ tru) (- rem divisor))
(values tru rem))))
(defknown %unary-ftruncate/single-float (single-float) single-float
(movable foldable flushable))
(defknown %unary-ftruncate/double-float (double-float) double-float
(movable foldable flushable))
(defknown %unary-ftruncate (real) float
(movable foldable flushable))
;; Convert (ftruncate x y) to the obvious implementation. We only
;; want this under certain conditions and let the generic ftruncate
;; handle the rest. (Note: if y = 1, the divide and multiply by y
;; should be removed by other deftransforms.)
(deftransform ftruncate ((x &optional (y 1))
(float &optional (or float integer)))
'(let ((res (%unary-ftruncate (/ x y))))
(values res (- x (* y res)))))
#+sparc
(defknown fast-unary-ftruncate ((or single-float double-float))
(or single-float double-float)
(movable foldable flushable))
#+sparc
(defoptimizer (fast-unary-ftruncate derive-type) ((f))
(one-arg-derive-type f
#'(lambda (n)
(ftruncate-derive-type-quot-aux n
(specifier-type '(integer 1 1))
nil))
#'ftruncate))
;; Convert %unary-ftruncate to unary-ftruncate/{single,double}-float
;; if x is known to be of the right type. Also, if the result is
;; known to fit in the same range as a (signed-byte 32), convert this
;; to %unary-truncate, which might be a single instruction, and float
;; the result. However, for sparc, we have a vop to do this so call
;; that, and for Sparc V9, we can actually handle a 64-bit integer
;; range.
(macrolet ((frob (ftype func)
`(deftransform %unary-ftruncate ((x) (,ftype))
(let* ((x-type (continuation-type x))
(lo (bound-value (numeric-type-low x-type)))
(hi (bound-value (numeric-type-high x-type)))
(limit-lo (- (ash 1 #-sparc-v9 31 #+sparc-v9 63)))
(limit-hi (ash 1 #-sparc-v9 31 #+sparc-v9 63)))
(if (and (numberp lo) (numberp hi)
(< limit-lo lo)
(< hi limit-hi))
#-sparc '(let ((result (coerce (%unary-truncate x) ',ftype)))
(if (zerop result)
(* result x)
result))
#+sparc '(let ((result (fast-unary-ftruncate x)))
(if (zerop result)
(* result x)
result))
'(,func x))))))
(frob single-float %unary-ftruncate/single-float)
(frob double-float %unary-ftruncate/double-float))
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
;;; FROUND
#-x87
(progn
(deftransform fround ((x &optional (y 1))
((or single-float double-float)
&optional (or single-float double-float integer)))
'(let ((res (%unary-fround (/ x y))))
(values res (- x (* y res)))))
(defknown %unary-fround (real) float
(movable foldable flushable))
(defknown %unary-fround/single-float (single-float) single-float
(movable foldable flushable))
(defknown %unary-fround/double-float (double-float) double-float
(movable foldable flushable))
(deftransform %unary-fround ((x) (single-float))
'(%unary-fround/single-float x))
(deftransform %unary-fround ((x) (double-float))
'(%unary-fround/double-float x))
); not x87
;;; Random:
;;;
(macrolet ((frob (fun type)
`(deftransform random ((num &optional state)
(,type &optional *) *
:when :both)
'(,fun num (or state *random-state*)))))
(frob %random-single-float single-float)
(frob %random-double-float double-float))
#-(or new-random random-mt19937)
(deftransform random ((num &optional state)
((integer 1 #.random-fixnum-max) &optional *))
_N"use inline fixnum operations"
;;; With the latest propagate-float-type code the compiler can inline
;;; truncate (signed-byte 32) allowing 31 bits, and (unsigned-byte 32)
;;; 32 bits on the x86. When not using the propagate-float-type
;;; feature the best size that can be inlined is 29 bits. The choice
;;; shouldn't cause bootstrap problems just slow code.
#+new-random
(deftransform random ((num &optional state)
((integer 1
&optional *))
#+x86 (intl:gettext "use inline (unsigned-byte 32) operations")
#-x86 (intl:gettext "use inline (signed-byte 32) operations")
'(values (truncate (%random-double-float (coerce num 'double-float)
(or state *random-state*)))))
#+random-mt19937
(deftransform random ((num &optional state)
((integer 1 #.(expt 2 32)) &optional *))
_N"use inline (unsigned-byte 32) operations"
(let* ((num-type (continuation-type num))
(num-high (cond ((numeric-type-p num-type)
(numeric-type-high num-type))
((union-type-p num-type)
;; Find the maximum of the union type. We
;; know this works because if we're in this
;; routine, NUM must be a subtype of
;; (INTEGER 1 2^32), so each member of the
;; union must be a subtype too.
(reduce #'max (union-type-types num-type)
:key #'numeric-type-high))
(t
(give-up)))))
;; Rather than doing (rem (random-chunk) num-high), we do,
;; essentially, (rem (* num-high (random-chunk)) #x100000000). I
;; (rtoy) believe this approach doesn't have the bias issue with
;; doing rem. This method works by treating (random-chunk) as if
;; it were a 32-bit fraction between 0 and 1, exclusive. Multiply
;; this by num-high to get a random number between 0 and num-high,
;; This should have no bias.
(cond ((constant-continuation-p num)
(if (= num-high (expt 2 32))
'(random-chunk (or state *random-state*))
'(values (bignum::%multiply
(random-chunk (or state *random-state*))
num))))
((< num-high (expt 2 32))
'(values (bignum::%multiply (random-chunk (or state *random-state*))
num)))
((= num-high (expt 2 32))
'(if (= num (expt 2 32))
(random-chunk (or state *random-state*))
(values (bignum::%multiply (random-chunk (or state *random-state*))
num))))
(error (intl:gettext "Shouldn't happen"))))))
;;;; Float accessors:
(defknown make-single-float ((signed-byte 32)) single-float
(movable foldable flushable))
(defknown make-double-float ((signed-byte 32) (unsigned-byte 32)) double-float
(movable foldable flushable))
(defknown single-float-bits (single-float) (signed-byte 32)
(movable foldable flushable))
(defknown double-float-high-bits (double-float) (signed-byte 32)
(movable foldable flushable))
(defknown double-float-low-bits (double-float) (unsigned-byte 32)
(movable foldable flushable))
(defknown double-float-bits (double-float)
(values (signed-byte 32) (unsigned-byte 32))
(movable foldable flushable))
#+double-double
(progn
(defknown double-double-float-p (t)
boolean
(movable foldable flushable))
(defknown %make-double-double-float (double-float double-float)
double-double-float
(movable foldable flushable))
(defknown double-double-hi (double-double-float)
double-float
(movable foldable flushable))
(defknown double-double-lo (double-double-float)
double-float
(movable foldable flushable))
) ; progn
(deftransform float-sign ((float &optional float2)
(single-float &optional single-float) *)
(if float2
(let ((temp (gensym)))
`(let ((,temp (abs float2)))
(if (minusp (single-float-bits float)) (- ,temp) ,temp)))
'(if (minusp (single-float-bits float)) -1f0 1f0)))
(deftransform float-sign ((float &optional float2)
(double-float &optional double-float) *)
(if float2
(let ((temp (gensym)))
`(let ((,temp (abs float2)))
(if (minusp (double-float-high-bits float)) (- ,temp) ,temp)))
'(if (minusp (double-float-high-bits float)) -1d0 1d0)))
#+double-double
(deftransform float-sign ((float &optional float2)
(double-double-float &optional double-double-float) *)
(if float2
(let ((temp (gensym)))
`(let ((,temp (abs float2)))
(if (minusp (float-sign (double-double-hi float)))
(- ,temp)
,temp)))
'(if (minusp (float-sign (double-double-hi float))) -1w0 1w0)))
;;;; DECODE-FLOAT, INTEGER-DECODE-FLOAT, SCALE-FLOAT:
;;;
;;; Convert these operations to format specific versions when the format is
(deftype single-float-exponent ()
`(integer ,(- vm:single-float-normal-exponent-min vm:single-float-bias
vm:single-float-digits)
,(- vm:single-float-normal-exponent-max vm:single-float-bias)))
(deftype double-float-exponent ()
`(integer ,(- vm:double-float-normal-exponent-min vm:double-float-bias
vm:double-float-digits)
,(- vm:double-float-normal-exponent-max vm:double-float-bias)))
(deftype single-float-int-exponent ()
`(integer ,(- vm:single-float-normal-exponent-min vm:single-float-bias
(* vm:single-float-digits 2))
,(- vm:single-float-normal-exponent-max vm:single-float-bias
vm:single-float-digits)))
`(integer ,(- vm:double-float-normal-exponent-min vm:double-float-bias
(* vm:double-float-digits 2))
,(- vm:double-float-normal-exponent-max vm:double-float-bias
vm:double-float-digits)))
`(integer 0 (,(ash 1 vm:single-float-digits))))
`(integer 0 (,(ash 1 vm:double-float-digits))))
(values (single-float 0.5f0 (1f0))
single-float-exponent
(member -1f0 1f0))
(movable foldable flushable))
(defknown decode-double-float (double-float)
(values (double-float 0.5d0 (1d0))
double-float-exponent
(member -1d0 1d0))
(movable foldable flushable))
(defknown integer-decode-single-float (single-float)
(values single-float-significand single-float-int-exponent (integer -1 1))
(movable foldable flushable))
(defknown integer-decode-double-float (double-float)
(values double-float-significand double-float-int-exponent (integer -1 1))
(defknown scale-single-float (single-float fixnum) single-float
(movable foldable flushable))
(defknown scale-double-float (double-float fixnum) double-float
(movable foldable flushable))
(deftransform decode-float ((x) (single-float) * :when :both)
(deftransform decode-float ((x) (double-float) * :when :both)
(deftransform integer-decode-float ((x) (single-float) * :when :both)
(deftransform integer-decode-float ((x) (double-float) * :when :both)
(deftransform scale-float ((f ex) (single-float *) * :when :both)
(if (and (backend-featurep :x86)
(not (backend-featurep :sse2))
(csubtypep (continuation-type ex)
(specifier-type '(signed-byte 32)))
(not (byte-compiling)))
'(coerce (%scalbn (coerce f 'double-float) ex) 'single-float)
'(scale-single-float f ex)))
(deftransform scale-float ((f ex) (double-float *) * :when :both)
(if (and (backend-featurep :x86)
(not (backend-featurep :sse2))
(csubtypep (continuation-type ex)
(specifier-type '(signed-byte 32))))
'(%scalbn f ex)
'(scale-double-float f ex)))
;;; toy@rtp.ericsson.se:
;;;
;;; Optimizers for scale-float. If the float has bounds, new bounds
;;; are computed for the result, if possible.
(defun scale-float-derive-type-aux (f ex same-arg)
(declare (ignore same-arg))
(flet ((scale-bound (x n)
;; We need to be a bit careful here and catch any overflows
;; that might occur. We can ignore underflows which become
;; zeros.
(set-bound
(let ((value (handler-case
(scale-float (bound-value x) n)
(floating-point-overflow ()
nil))))
;; This check is necessary for ppc because the current
;; implementation on ppc doesn't signal floating-point
;; overflow. (How many other places do we need to check
;; for this?)
(if (and (floatp value) (float-infinity-p value))
nil
value))
(consp x))))
(when (and (numeric-type-p f) (numeric-type-p ex))
(let ((f-lo (numeric-type-low f))
(f-hi (numeric-type-high f))
(ex-lo (numeric-type-low ex))
(ex-hi (numeric-type-high ex))
(new-lo nil)
(new-hi nil))
(when (and f-hi ex-hi)
(setf new-hi (scale-bound f-hi ex-hi)))
(when (and f-lo ex-lo)
(setf new-lo (scale-bound f-lo ex-lo)))
;; We're computing bounds for scale-float. Assume the bounds
;; on f are fl and fh, and the bounds on ex are nl and nh.
;; The resulting bound should be fl*2^nl and fh*2^nh.
;; However, if fh is negative, and we get an underflow, we
;; might get bounds like 0 and fh*2^nh < 0. Our bounds are
;; backwards. Thus, swap the bounds to get the correct
;; bounds.
(when (and new-lo new-hi (< (bound-value new-hi)
(bound-value new-lo)))
(rotatef new-lo new-hi))
(make-numeric-type :class (numeric-type-class f)
:format (numeric-type-format f)
:complexp :real
:low new-lo
:high new-hi)))))
;;;
(defoptimizer (scale-float derive-type) ((f ex))
(two-arg-derive-type f ex #'scale-float-derive-type-aux
#'scale-float t))
;;; Defoptimizers for %single-float and %double-float. This makes the
;;; FLOAT function return the correct ranges if the input has some
;;; defined range. Quite useful if we want to convert some type of
(macrolet
((frob (fun type)
(let ((aux-name (symbolicate fun "-DERIVE-TYPE-AUX")))
(defun ,aux-name (num)
;; When converting a number to a float, the limits are
;; the "same."
(let* ((lo (bound-func #'(lambda (x)
;; If we can't coerce it, we
;; return a NIL for the bound.
;; (Is IGNORE-ERRORS too
;; heavy-handed? Should we
;; try to do something more
;; fine-grained?)
(ignore-errors (coerce x ',type)))
(numeric-type-low num)))
(hi (bound-func #'(lambda (x)
(ignore-errors (coerce x ',type)))
(numeric-type-high num))))
(specifier-type `(,',type ,(or lo '*) ,(or hi '*)))))
(defoptimizer (,fun derive-type) ((num))
(frob %single-float single-float)
(frob %double-float double-float))
;;;; Float contagion:
;;; FLOAT-CONTAGION-ARG1, ARG2 -- Internal
;;;
;;; Do some stuff to recognize when the loser is doing mixed float and
;;; rational arithmetic, or different float types, and fix it up. If we don't,
;;; he won't even get so much as an efficency note.
;;;
(deftransform float-contagion-arg1 ((x y) * * :defun-only t :node node)
`(,(continuation-function-name (basic-combination-fun node))
(float x y) y))
;;;
(deftransform float-contagion-arg2 ((x y) * * :defun-only t :node node)
`(,(continuation-function-name (basic-combination-fun node))
x (float y x)))
(dolist (x '(+ * / -))
(%deftransform x '(function (rational float) *) #'float-contagion-arg1)
(%deftransform x '(function (float rational) *) #'float-contagion-arg2))
(dolist (x '(= < > + * / -))
(%deftransform x '(function (single-float double-float) *)
#'float-contagion-arg1)
(%deftransform x '(function (double-float single-float) *)
#'float-contagion-arg2))
;;; Prevent zerop, plusp, minusp from losing horribly. We can't in general
;;; float rational args to comparison, since Common Lisp semantics says we are
;;; supposed to compare as rationals, but we can do it for any rational that
;;; has a precise representation as a float (such as 0).
`(deftransform ,op ((x y) (float rational) * :when :both)
(unless (constant-continuation-p y)
(give-up (intl:gettext "Can't open-code float to rational comparison.")))
(let ((val (continuation-value y)))
(unless (eql (rational (float val)) val)
(give-up (intl:gettext "~S doesn't have a precise float representation.")
val)))
`(,',op x (float y x)))))
;; Convert (/ x n) to (* x (/ n)) when x is a float and n is a power
;; of two, because (/ n) can be reprsented exactly.
(deftransform / ((x y) (float float) * :when :both)
(unless (constant-continuation-p y)
(give-up))
(let ((val (continuation-value y)))
(multiple-value-bind (frac exp sign)
(decode-float val)
(unless (= frac 0.5)
(give-up))
`(* x (float (/ ,val) x)))))
;; Convert 2*x to x+x.
(deftransform * ((x y) (float real) * :when :both)
(unless (constant-continuation-p y)
(give-up))
(let ((val (continuation-value y)))
(unless (= val 2)
(give-up))
'(+ x x)))
;;;; Irrational transforms:
(defknown (%tan %sinh %asinh %atanh %log %logb %log10 #+x87 %tan-quick)
(double-float) double-float
(movable foldable flushable))
(defknown (%sin %cos %tanh #+x87 %sin-quick #+x87 %cos-quick)
(double-float) (double-float -1.0d0 1.0d0)
(movable foldable flushable))
(defknown (%asin %atan)
(double-float) (double-float #.(- (/ pi 2)) #.(/ pi 2))
(movable foldable flushable))
(defknown (%acos)
(double-float) (double-float 0.0d0 #.pi)
(movable foldable flushable))
(defknown (%cosh)
(double-float) (double-float 1.0d0)
(movable foldable flushable))
(defknown (%acosh %exp %sqrt)
(double-float) (double-float 0.0d0)
(movable foldable flushable))
(defknown %expm1
(double-float) (double-float -1d0)
(movable foldable flushable))
(defknown (%hypot)
(double-float double-float) (double-float 0d0)
(movable foldable flushable))
(defknown (%pow)
(double-float double-float) double-float
(movable foldable flushable))
(defknown (%atan2)
(double-float double-float) (double-float #.(- pi) #.pi)
(movable foldable flushable))
(defknown (%scalb)
(double-float double-float) double-float
(movable foldable flushable))
(defknown (%scalbn)
(double-float (signed-byte 32)) double-float
(movable foldable flushable))
(double-float) double-float
(movable foldable flushable))
(dolist (stuff '((exp %exp *)
(log %log float)
(sqrt %sqrt float)
(asin %asin float)
(acos %acos float)
(atan %atan *)
(sinh %sinh *)
(cosh %cosh *)
(tanh %tanh *)
(asinh %asinh *)
(acosh %acosh float)
(atanh %atanh float)))
(destructuring-bind (name prim rtype) stuff
(deftransform name ((x) '(single-float) rtype :eval-name t)
`(coerce (,prim (coerce x 'double-float)) 'single-float))
(deftransform name ((x) '(double-float) rtype :eval-name t :when :both)
`(,prim x))))
;;; The argument range is limited on the x86 FP trig. functions. A
;;; post-test can detect a failure (and load a suitable result), but
;;; this test is avoided if possible.
;;;
;;; Simple tests show that sin/cos produce numbers greater than 1 when
;;; the arg >= 2^63. tan produces floating-point invalid exceptions
;;; for arg >= 2^62. So limit these to that range.
(dolist (stuff '((sin %sin %sin-quick 63)
(cos %cos %cos-quick 63)
(tan %tan %tan-quick 62)))
(destructuring-bind (name prim prim-quick limit)
stuff
(deftransform name ((x) '(single-float) '* :eval-name t)
(if (and (backend-featurep :x86)
(not (backend-featurep :sse2)))
(cond ((csubtypep (continuation-type x)
(specifier-type `(single-float
(,(- (expt 2f0 limit)))
(,(expt 2f0 limit)))))
`(coerce (,prim-quick (coerce x 'double-float))
'single-float))
(t
(compiler-note
_N"Unable to avoid inline argument range check~@
because the argument range (~s) was not within 2^~D"
(type-specifier (continuation-type x))
limit)
`(coerce (,prim (coerce x 'double-float)) 'single-float)))
`(coerce (,prim (coerce x 'double-float)) 'single-float)))
(deftransform name ((x) '(double-float) '* :eval-name t :when :both)
(if (and (backend-featurep :x86)
(not (backend-featurep :sse2)))
(cond ((csubtypep (continuation-type x)
(specifier-type `(double-float
(,(- (expt 2d0 limit)))
(,(expt 2d0 limit)))))
`(,prim-quick x))
(t
(compiler-note
_N"Unable to avoid inline argument range check~@
because the argument range (~s) was not within 2^~D"
(type-specifier (continuation-type x))
limit)
`(,prim x)))
`(,prim x)))))
(deftransform atan ((x y) (single-float single-float) *)
`(coerce (%atan2 (coerce x 'double-float) (coerce y 'double-float))
'single-float))
(deftransform atan ((x y) (double-float double-float) * :when :both)
`(%atan2 x y))
(deftransform expt ((x y) ((single-float 0f0) single-float) *)
`(coerce (%pow (coerce x 'double-float) (coerce y 'double-float))
'single-float))
(deftransform expt ((x y) ((double-float 0d0) double-float) * :when :both)
`(%pow x y))
(deftransform expt ((x y) ((single-float 0f0) (signed-byte 32)) *)
`(coerce (%pow (coerce x 'double-float) (coerce y 'double-float))
'single-float))
(deftransform expt ((x y) ((double-float 0d0) (signed-byte 32)) * :when :both)
`(%pow x (coerce y 'double-float)))
;;; ANSI says log with base zero returns zero.
(deftransform log ((x y) (float float) float)
'(if (zerop y) y (/ (log x) (log y))))
;;; Handle some simple transformations
(deftransform abs ((x) ((complex double-float)) double-float :when :both)
'(%hypot (realpart x) (imagpart x)))
(deftransform abs ((x) ((complex single-float)) single-float)
'(coerce (%hypot (coerce (realpart x) 'double-float)
(coerce (imagpart x) 'double-float))
(deftransform abs ((x) (real) real)
(let ((x-type (continuation-type x)))
;; If the arg is known to non-negative, we can just return the
;; arg. However, (abs -0.0) is 0.0, so this transform only works
;; on floats that are known not to include negative zero.
(if (csubtypep x-type (specifier-type '(or (rational 0) (float (0d0)) (member 0f0 0d0))))
'x
(give-up))))
(deftransform phase ((x) ((complex double-float)) double-float :when :both)
'(%atan2 (imagpart x) (realpart x)))
(deftransform phase ((x) ((complex single-float)) single-float)
'(coerce (%atan2 (coerce (imagpart x) 'double-float)
(coerce (realpart x) 'double-float))
'single-float))
(deftransform phase ((x) ((float)) float :when :both)
'(if (minusp (float-sign x))
;;; The number is of type REAL.
(defun numeric-type-real-p (type)
(and (numeric-type-p type)
(eq (numeric-type-complexp type) :real)))
;;; Coerce a numeric type bound to the given type while handling
;;; exclusive bounds.
(defun coerce-numeric-bound (bound type)
(when bound
(if (consp bound)
(list (coerce (car bound) type))
(coerce bound type))))
;;;; Optimizers for elementary functions
;;;;
;;;; These optimizers compute the output range of the elementary
;;;; function, based on the domain of the input.
;;;;
;;; Generate a specifier for a complex type specialized to the same
;;; type as the argument.
(defun complex-float-type (arg)
(declare (type numeric-type arg))
(let* ((format (case (numeric-type-class arg)
((integer rational) 'single-float)
(t (numeric-type-format arg))))
(float-type (or format 'float)))
(specifier-type `(complex ,float-type))))
;;; Compute a specifier like '(or float (complex float)), except float
;;; should be the right kind of float. Allow bounds for the float
;;; part too.
(defun float-or-complex-float-type (arg &optional lo hi)
(declare (type numeric-type arg))
(let* ((format (case (numeric-type-class arg)
((integer rational) 'single-float)
(t (numeric-type-format arg))))
(float-type (or format 'float))
(lo (coerce-numeric-bound lo float-type))
(hi (coerce-numeric-bound hi float-type)))
(specifier-type `(or (,float-type ,(or lo '*) ,(or hi '*))
(complex ,float-type)))))
;;; Domain-Subtype
;;;
;;; Test if the numeric-type ARG is within in domain specified by
;;; DOMAIN-LOW and DOMAIN-HIGH, consider negative and positive zero to
;;; be distinct as for the :negative-zero-is-not-zero feature. Note
;;; that only inclusive and open domain limits are handled as these
;;; are the only types of limits currently used. With the
;;; :negative-zero-is-not-zero feature this could be handled by the
;;; numeric subtype code in type.lisp.
;;;
(defun domain-subtypep (arg domain-low domain-high)
(declare (type numeric-type arg)
(type (or real null) domain-low domain-high))
(let* ((arg-lo (numeric-type-low arg))
(arg-lo-val (bound-value arg-lo))
(arg-hi (numeric-type-high arg))
(arg-hi-val (bound-value arg-hi)))
;; Check that the ARG bounds are correctly canonicalised.
(when (and arg-lo (floatp arg-lo-val) (zerop arg-lo-val) (consp arg-lo)
(minusp (float-sign arg-lo-val)))
(compiler-note _N"Float zero bound ~s not correctly canonicalised?" arg-lo)
(when (and arg-hi (zerop arg-hi-val) (floatp arg-hi-val) (consp arg-hi)
(plusp (float-sign arg-hi-val)))
(compiler-note _N"Float zero bound ~s not correctly canonicalised?" arg-hi)
(setq arg-hi -0l0 arg-hi-val -0l0))
(flet ((fp-neg-zero-p (f) ; Is F -0.0?
(and (floatp f) (zerop f) (minusp (float-sign f))))
(fp-pos-zero-p (f) ; Is F +0.0?
(and (floatp f) (zerop f) (plusp (float-sign f)))))
(and (or (null domain-low)
(and arg-lo (>= arg-lo-val domain-low)
(not (and (fp-pos-zero-p domain-low)
(fp-neg-zero-p arg-lo)))))
(or (null domain-high)
(and arg-hi (<= arg-hi-val domain-high)
(not (and (fp-neg-zero-p domain-high)
(fp-pos-zero-p arg-hi)))))))))
;;; Handle monotonic functions of a single variable whose domain is
;;; possibly part of the real line. ARG is the variable, FCN is the
;;; function, and DOMAIN is a specifier that gives the (real) domain
;;; of the function. If ARG is a subset of the DOMAIN, we compute the
;;; bounds directly. Otherwise, we compute the bounds for the
;;; intersection between ARG and DOMAIN, and then append a complex
;;; result, which occurs for the parts of ARG not in the DOMAIN.
;;;
;;; Negative and positive zero are considered distinct within
;;; DOMAIN-LOW and DOMAIN-HIGH, as for the :negative-zero-is-not-zero
;;; feature.
;;;
;;; DEFAULT-LOW and DEFAULT-HIGH are the lower and upper bounds if we
;;; can't compute the bounds using FCN.
;;;
(defun elfun-derive-type-simple (arg fcn domain-low domain-high
default-low default-high
(declare (type (or null real) domain-low domain-high))
(etypecase arg
(numeric-type
(cond ((eq (numeric-type-complexp arg) :complex)
(complex-float-type arg))
((numeric-type-real-p arg)
;; The argument is real, so let's find the intersection
;; between the argument and the domain of the function.
;; We compute the bounds on the intersection, and for
;; everything else, we return a complex number of the
;; appropriate type.
(multiple-value-bind (intersection difference)
(interval-intersection/difference
(numeric-type->interval arg)
(make-interval :low domain-low :high domain-high))
(cond
(intersection
;; Process the intersection.
(let* ((low (interval-low intersection))
(high (interval-high intersection))
(res-lo (or (bound-func fcn (if increasingp low high))
(res-hi (or (bound-func fcn (if increasingp high low))
;; Result specifier type.
(format (case (numeric-type-class arg)
((integer rational) 'single-float)
(t (numeric-type-format arg))))
(bound-type (or format 'float))
(result-type
(make-numeric-type
:class 'float
:format format
:low (coerce-numeric-bound res-lo bound-type)
:high (coerce-numeric-bound res-hi bound-type))))
;; If the ARG is a subset of the domain, we don't
;; have to worry about the difference, because that
;; can't occur.
(if (or (null difference)
;; Check if the arg is within the domain.
(domain-subtypep arg domain-low domain-high))
result-type
(list result-type
(specifier-type `(complex ,bound-type))))))
(t
;; No intersection so the result must be purely complex.
(complex-float-type arg)))))
(float-or-complex-float-type arg default-low default-high))))))
((frob (name domain-low domain-high def-low-bnd def-high-bnd
&key (increasingp t))
(let ((num (gensym)))
`(defoptimizer (,name derive-type) ((,num))
(one-arg-derive-type
,num
#'(lambda (arg)
(elfun-derive-type-simple arg #',name
,domain-low ,domain-high
,def-low-bnd ,def-high-bnd
;; These functions are easy because they are defined for the whole