X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcompiler%2Ffloat-tran.lisp;h=882cc7004e298e9f28709cabdb8bd0ee989a0f79;hb=d335afdcf50b641a0aafd32741e0777d4e12a59b;hp=262f86777cae2a22f469e98076068ceda5b4af71;hpb=04f74c8cbe98704aa761e187741984bc8fe3983f;p=sbcl.git diff --git a/src/compiler/float-tran.lisp b/src/compiler/float-tran.lisp index 262f867..882cc70 100644 --- a/src/compiler/float-tran.lisp +++ b/src/compiler/float-tran.lisp @@ -15,8 +15,10 @@ ;;;; coercions -(defknown %single-float (real) single-float (movable foldable)) -(defknown %double-float (real) double-float (movable foldable)) +(defknown %single-float (real) single-float + (movable foldable)) +(defknown %double-float (real) double-float + (movable foldable)) (deftransform float ((n f) (* single-float) *) '(%single-float n)) @@ -332,6 +334,73 @@ (%deftransform x '(function (double-float single-float) *) #'float-contagion-arg2)) +(macrolet ((def (type &rest args) + `(deftransform * ((x y) (,type (constant-arg (member ,@args))) * + ;; Beware the SNaN! + :policy (zerop float-accuracy)) + "optimize multiplication by one" + (let ((y (lvar-value y))) + (if (minusp y) + '(%negate x) + 'x))))) + (def single-float 1.0 -1.0) + (def double-float 1.0d0 -1.0d0)) + +;;; Return the reciprocal of X if it can be represented exactly, NIL otherwise. +(defun maybe-exact-reciprocal (x) + (unless (zerop x) + (handler-case + (multiple-value-bind (significand exponent sign) + (integer-decode-float x) + ;; only powers of 2 can be inverted exactly + (unless (zerop (logand significand (1- significand))) + (return-from maybe-exact-reciprocal nil)) + (let ((expected (/ sign significand (expt 2 exponent))) + (reciprocal (/ x))) + (multiple-value-bind (significand exponent sign) + (integer-decode-float reciprocal) + ;; Denorms can't be inverted safely. + (and (eql expected (* sign significand (expt 2 exponent))) + reciprocal)))) + (error () (return-from maybe-exact-reciprocal nil))))) + +;;; Replace constant division by multiplication with exact reciprocal, +;;; if one exists. +(macrolet ((def (type) + `(deftransform / ((x y) (,type (constant-arg ,type)) * + :node node) + "convert to multiplication by reciprocal" + (let ((n (lvar-value y))) + (if (policy node (zerop float-accuracy)) + `(* x ,(/ n)) + (let ((r (maybe-exact-reciprocal n))) + (if r + `(* x ,r) + (give-up-ir1-transform + "~S does not have an exact reciprocal" + n)))))))) + (def single-float) + (def double-float)) + +;;; Optimize addition and subtraction of zero +(macrolet ((def (op type &rest args) + `(deftransform ,op ((x y) (,type (constant-arg (member ,@args))) * + ;; Beware the SNaN! + :policy (zerop float-accuracy)) + 'x))) + ;; No signed zeros, thanks. + (def + single-float 0 0.0) + (def - single-float 0 0.0) + (def + double-float 0 0.0 0.0d0) + (def - double-float 0 0.0 0.0d0)) + +;;; On most platforms (+ x x) is faster than (* x 2) +(macrolet ((def (type &rest args) + `(deftransform * ((x y) (,type (constant-arg (member ,@args)))) + '(+ x x)))) + (def single-float 2 2.0) + (def double-float 2 2.0 2.0d0)) + ;;; Prevent ZEROP, PLUSP, and 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 @@ -468,27 +537,27 @@ (deftransform ,name ((x) (single-float) *) #!+x86 (cond ((csubtypep (lvar-type x) (specifier-type '(single-float - (#.(- (expt 2f0 64))) - (#.(expt 2f0 64))))) + (#.(- (expt 2f0 63))) + (#.(expt 2f0 63))))) `(coerce (,',prim-quick (coerce x 'double-float)) 'single-float)) (t (compiler-notify "unable to avoid inline argument range check~@ - because the argument range (~S) was not within 2^64" + because the argument range (~S) was not within 2^63" (type-specifier (lvar-type x))) `(coerce (,',prim (coerce x 'double-float)) 'single-float))) #!-x86 `(coerce (,',prim (coerce x 'double-float)) 'single-float)) (deftransform ,name ((x) (double-float) *) #!+x86 (cond ((csubtypep (lvar-type x) (specifier-type '(double-float - (#.(- (expt 2d0 64))) - (#.(expt 2d0 64))))) + (#.(- (expt 2d0 63))) + (#.(expt 2d0 63))))) `(,',prim-quick x)) (t (compiler-notify "unable to avoid inline argument range check~@ - because the argument range (~S) was not within 2^64" + because the argument range (~S) was not within 2^63" (type-specifier (lvar-type x))) `(,',prim x))) #!-x86 `(,',prim x))))) @@ -636,7 +705,7 @@ (progn ;;; Handle monotonic functions of a single variable whose domain is -;;; possibly part of the real line. ARG is the variable, FCN is the +;;; possibly part of the real line. ARG is the variable, FUN 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 @@ -647,8 +716,8 @@ ;;; DOMAIN-LOW and DOMAIN-HIGH. ;;; ;;; 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 +;;; can't compute the bounds using FUN. +(defun elfun-derive-type-simple (arg fun domain-low domain-high default-low default-high &optional (increasingp t)) (declare (type (or null real) domain-low domain-high)) @@ -672,9 +741,9 @@ ;; Process the intersection. (let* ((low (interval-low intersection)) (high (interval-high intersection)) - (res-lo (or (bound-func fcn (if increasingp low high)) + (res-lo (or (bound-func fun (if increasingp low high)) default-low)) - (res-hi (or (bound-func fcn (if increasingp high low)) + (res-hi (or (bound-func fun (if increasingp high low)) default-high)) (format (case (numeric-type-class arg) ((integer rational) 'single-float) @@ -953,7 +1022,7 @@ ;; But a positive real to any power is well-defined. (merged-interval-expt x y)) ((and (csubtypep x (specifier-type 'rational)) - (csubtypep x (specifier-type 'rational))) + (csubtypep y (specifier-type 'rational))) ;; A rational to the power of a rational could be a rational ;; or a possibly-complex single float (specifier-type '(or rational single-float (complex single-float)))) @@ -1178,27 +1247,45 @@ ;;; of complex operation VOPs. (macrolet ((frob (type) `(progn + (deftransform complex ((r) (,type)) + '(complex r ,(coerce 0 type))) + (deftransform complex ((r i) (,type (and real (not ,type)))) + '(complex r (truly-the ,type (coerce i ',type)))) + (deftransform complex ((r i) ((and real (not ,type)) ,type)) + '(complex (truly-the ,type (coerce r ',type)) i)) ;; negation + #!-complex-float-vops (deftransform %negate ((z) ((complex ,type)) *) '(complex (%negate (realpart z)) (%negate (imagpart z)))) ;; complex addition and subtraction + #!-complex-float-vops (deftransform + ((w z) ((complex ,type) (complex ,type)) *) '(complex (+ (realpart w) (realpart z)) (+ (imagpart w) (imagpart z)))) + #!-complex-float-vops (deftransform - ((w z) ((complex ,type) (complex ,type)) *) '(complex (- (realpart w) (realpart z)) (- (imagpart w) (imagpart z)))) ;; Add and subtract a complex and a real. + #!-complex-float-vops (deftransform + ((w z) ((complex ,type) real) *) - '(complex (+ (realpart w) z) (imagpart w))) + `(complex (+ (realpart w) z) + (+ (imagpart w) ,(coerce 0 ',type)))) + #!-complex-float-vops (deftransform + ((z w) (real (complex ,type)) *) - '(complex (+ (realpart w) z) (imagpart w))) + `(complex (+ (realpart w) z) + (+ (imagpart w) ,(coerce 0 ',type)))) ;; Add and subtract a real and a complex number. + #!-complex-float-vops (deftransform - ((w z) ((complex ,type) real) *) - '(complex (- (realpart w) z) (imagpart w))) + `(complex (- (realpart w) z) + (- (imagpart w) ,(coerce 0 ',type)))) + #!-complex-float-vops (deftransform - ((z w) (real (complex ,type)) *) - '(complex (- z (realpart w)) (- (imagpart w)))) + `(complex (- z (realpart w)) + (- ,(coerce 0 ',type) (imagpart w)))) ;; Multiply and divide two complex numbers. + #!-complex-float-vops (deftransform * ((x y) ((complex ,type) (complex ,type)) *) '(let* ((rx (realpart x)) (ix (imagpart x)) @@ -1207,39 +1294,81 @@ (complex (- (* rx ry) (* ix iy)) (+ (* rx iy) (* ix ry))))) (deftransform / ((x y) ((complex ,type) (complex ,type)) *) + #!-complex-float-vops '(let* ((rx (realpart x)) (ix (imagpart x)) (ry (realpart y)) (iy (imagpart y))) (if (> (abs ry) (abs iy)) (let* ((r (/ iy ry)) - (dn (* ry (+ 1 (* r r))))) + (dn (+ ry (* r iy)))) (complex (/ (+ rx (* ix r)) dn) (/ (- ix (* rx r)) dn))) (let* ((r (/ ry iy)) - (dn (* iy (+ 1 (* r r))))) + (dn (+ iy (* r ry)))) (complex (/ (+ (* rx r) ix) dn) - (/ (- (* ix r) rx) dn)))))) + (/ (- (* ix r) rx) dn))))) + #!+complex-float-vops + `(let* ((cs (conjugate (sb!vm::swap-complex x))) + (ry (realpart y)) + (iy (imagpart y))) + (if (> (abs ry) (abs iy)) + (let* ((r (/ iy ry)) + (dn (+ ry (* r iy)))) + (/ (+ x (* cs r)) dn)) + (let* ((r (/ ry iy)) + (dn (+ iy (* r ry)))) + (/ (+ (* x r) cs) dn))))) ;; Multiply a complex by a real or vice versa. + #!-complex-float-vops (deftransform * ((w z) ((complex ,type) real) *) '(complex (* (realpart w) z) (* (imagpart w) z))) + #!-complex-float-vops (deftransform * ((z w) (real (complex ,type)) *) '(complex (* (realpart w) z) (* (imagpart w) z))) - ;; Divide a complex by a real. + ;; Divide a complex by a real or vice versa. + #!-complex-float-vops (deftransform / ((w z) ((complex ,type) real) *) '(complex (/ (realpart w) z) (/ (imagpart w) z))) + (deftransform / ((x y) (,type (complex ,type)) *) + #!-complex-float-vops + '(let* ((ry (realpart y)) + (iy (imagpart y))) + (if (> (abs ry) (abs iy)) + (let* ((r (/ iy ry)) + (dn (+ ry (* r iy)))) + (complex (/ x dn) + (/ (- (* x r)) dn))) + (let* ((r (/ ry iy)) + (dn (+ iy (* r ry)))) + (complex (/ (* x r) dn) + (/ (- x) dn))))) + #!+complex-float-vops + '(let* ((ry (realpart y)) + (iy (imagpart y))) + (if (> (abs ry) (abs iy)) + (let* ((r (/ iy ry)) + (dn (+ ry (* r iy)))) + (/ (complex x (- (* x r))) dn)) + (let* ((r (/ ry iy)) + (dn (+ iy (* r ry)))) + (/ (complex (* x r) (- x)) dn))))) ;; conjugate of complex number + #!-complex-float-vops (deftransform conjugate ((z) ((complex ,type)) *) '(complex (realpart z) (- (imagpart z)))) ;; CIS (deftransform cis ((z) ((,type)) *) '(complex (cos z) (sin z))) ;; comparison + #!-complex-float-vops (deftransform = ((w z) ((complex ,type) (complex ,type)) *) '(and (= (realpart w) (realpart z)) (= (imagpart w) (imagpart z)))) + #!-complex-float-vops (deftransform = ((w z) ((complex ,type) real) *) '(and (= (realpart w) z) (zerop (imagpart w)))) + #!-complex-float-vops (deftransform = ((w z) (real (complex ,type)) *) '(and (= (realpart z) w) (zerop (imagpart z))))))) @@ -1253,7 +1382,7 @@ ;;; inputs are union types. #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.) (progn -(defun trig-derive-type-aux (arg domain fcn +(defun trig-derive-type-aux (arg domain fun &optional def-lo def-hi (increasingp t)) (etypecase arg (numeric-type @@ -1274,8 +1403,8 @@ ;; exactly the same way as the functions themselves do ;; it. (if (csubtypep arg domain) - (let ((res-lo (bound-func fcn (numeric-type-low arg))) - (res-hi (bound-func fcn (numeric-type-high arg)))) + (let ((res-lo (bound-func fun (numeric-type-low arg))) + (res-hi (bound-func fun (numeric-type-high arg)))) (unless increasingp (rotatef res-lo res-hi)) (make-numeric-type @@ -1370,15 +1499,51 @@ (define-frobs truncate %unary-truncate) (define-frobs round %unary-round)) -;;; 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 ((defaulted-y (if y 'y 1))) - `(let ((res (%unary-truncate (/ x ,defaulted-y)))) - (values res (- x (* ,defaulted-y res)))))) +(deftransform %unary-truncate ((x) (single-float)) + `(%unary-truncate/single-float x)) +(deftransform %unary-truncate ((x) (double-float)) + `(%unary-truncate/double-float x)) + +;;; Convert (TRUNCATE x y) to the obvious implementation. +;;; +;;; ...plus hair: Insert explicit coercions to appropriate float types: Python +;;; is reluctant it generate explicit integer->float coercions due to +;;; precision issues (see SAFE-SINGLE-COERCION-P &co), but this is not an +;;; issue here as there is no DERIVE-TYPE optimizer on specialized versions of +;;; %UNARY-TRUNCATE, so the derived type of TRUNCATE remains the best we can +;;; do here -- which is fine. Also take care not to add unnecassary division +;;; or multiplication by 1, since we are not able to always eliminate them, +;;; depending on FLOAT-ACCURACY. Finally, leave out the secondary value when +;;; we know it is unused: COERCE is not flushable. +(macrolet ((def (type other-float-arg-types) + (let ((unary (symbolicate "%UNARY-TRUNCATE/" type)) + (coerce (symbolicate "%" type))) + `(deftransform truncate ((x &optional y) + (,type + &optional (or ,type ,@other-float-arg-types integer)) + * :result result) + (let* ((result-type (and result + (lvar-derived-type result))) + (compute-all (and (values-type-p result-type) + (not (type-single-value-p result-type))))) + (if (or (not y) + (and (constant-lvar-p y) (= 1 (lvar-value y)))) + (if compute-all + `(let ((res (,',unary x))) + (values res (- x (,',coerce res)))) + `(let ((res (,',unary x))) + ;; Dummy secondary value! + (values res x))) + (if compute-all + `(let* ((f (,',coerce y)) + (res (,',unary (/ x f)))) + (values res (- x (* f (,',coerce res))))) + `(let* ((f (,',coerce y)) + (res (,',unary (/ x f)))) + ;; Dummy secondary value! + (values res x))))))))) + (def single-float ()) + (def double-float (single-float))) (deftransform floor ((number &optional divisor) (float &optional (or integer float)))