X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcompiler%2Ffloat-tran.lisp;h=a1184e12ae59221f26aa7c22738b823fe6ca8253;hb=95591ed483dbb8c0846c129953acac1554f28809;hp=1cb39f6654687cee0158bda7ea1680344fa45501;hpb=0a7604d54581d2c846838c26ce6a7993629586fa;p=sbcl.git

diff --git a/src/compiler/float-tran.lisp b/src/compiler/float-tran.lisp
index 1cb39f6..a1184e1 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 flushable))
-(defknown %double-float (real) double-float (movable foldable flushable))
+(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))
@@ -43,6 +45,59 @@
                 '(,fun num (or state *random-state*)))))
   (frob %random-single-float single-float)
   (frob %random-double-float double-float))
+
+;;; Mersenne Twister RNG
+;;;
+;;; FIXME: It's unpleasant to have RANDOM functionality scattered
+;;; through the code this way. It would be nice to move this into the
+;;; same file as the other RANDOM definitions.
+(deftransform random ((num &optional state)
+                      ((integer 1 #.(expt 2 sb!vm::n-word-bits)) &optional *))
+  ;; FIXME: I almost conditionalized this as #!+sb-doc. Find some way
+  ;; of automatically finding #!+sb-doc in proximity to DEFTRANSFORM
+  ;; to let me scan for places that I made this mistake and didn't
+  ;; catch myself.
+  "use inline (UNSIGNED-BYTE 32) operations"
+  (let ((type (lvar-type num))
+        (limit (expt 2 sb!vm::n-word-bits))
+        (random-chunk (ecase sb!vm::n-word-bits
+                        (32 'random-chunk)
+                        (64 'sb!kernel::big-random-chunk))))
+    (if (numeric-type-p type)
+        (let ((num-high (numeric-type-high (lvar-type num))))
+          (aver num-high)
+          (cond ((constant-lvar-p num)
+                 ;; Check the worst case sum absolute error for the
+                 ;; random number expectations.
+                 (let ((rem (rem limit num-high)))
+                   (unless (< (/ (* 2 rem (- num-high rem))
+                                 num-high limit)
+                              (expt 2 (- sb!kernel::random-integer-extra-bits)))
+                     (give-up-ir1-transform
+                      "The random number expectations are inaccurate."))
+                   (if (= num-high limit)
+                       `(,random-chunk (or state *random-state*))
+                       #!-(or x86 x86-64)
+                       `(rem (,random-chunk (or state *random-state*)) num)
+                       #!+(or x86 x86-64)
+                       ;; Use multiplication, which is faster.
+                       `(values (sb!bignum::%multiply
+                                 (,random-chunk (or state *random-state*))
+                                 num)))))
+                ((> num-high random-fixnum-max)
+                 (give-up-ir1-transform
+                  "The range is too large to ensure an accurate result."))
+                #!+(or x86 x86-64)
+                ((< num-high limit)
+                 `(values (sb!bignum::%multiply
+                           (,random-chunk (or state *random-state*))
+                           num)))
+                (t
+                 `(rem (,random-chunk (or state *random-state*)) num))))
+        ;; KLUDGE: a relatively conservative treatment, but better
+        ;; than a bug (reported by PFD sbcl-devel towards the end of
+        ;; 2004-11.
+        '(rem (random-chunk (or state *random-state*)) num))))
 
 ;;;; float accessors
 
@@ -233,7 +288,10 @@
                 (specifier-type `(,',type ,(or lo '*) ,(or hi '*)))))
 
             (defoptimizer (,fun derive-type) ((num))
-              (one-arg-derive-type num #',aux-name #',fun))))))
+              (handler-case
+                  (one-arg-derive-type num #',aux-name #',fun)
+                (type-error ()
+                  nil)))))))
   (frob %single-float single-float
         most-negative-single-float most-positive-single-float)
   (frob %double-float double-float
@@ -242,15 +300,29 @@
 
 ;;;; float contagion
 
+(defun safe-ctype-for-single-coercion-p (x)
+  ;; See comment in SAFE-SINGLE-COERCION-P -- this deals with the same
+  ;; problem, but in the context of evaluated and compiled (+ <int> <single>)
+  ;; giving different result if we fail to check for this.
+  (or (not (csubtypep x (specifier-type 'integer)))
+      (csubtypep x (specifier-type `(integer ,most-negative-exactly-single-float-fixnum
+                                             ,most-positive-exactly-single-float-fixnum)))))
+
 ;;; 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 efficiency note.
 (deftransform float-contagion-arg1 ((x y) * * :defun-only t :node node)
-  `(,(lvar-fun-name (basic-combination-fun node))
-    (float x y) y))
+  (if (or (not (types-equal-or-intersect (lvar-type y) (specifier-type 'single-float)))
+          (safe-ctype-for-single-coercion-p (lvar-type x)))
+      `(,(lvar-fun-name (basic-combination-fun node))
+         (float x y) y)
+      (give-up-ir1-transform)))
 (deftransform float-contagion-arg2 ((x y) * * :defun-only t :node node)
-  `(,(lvar-fun-name (basic-combination-fun node))
-    x (float y x)))
+  (if (or (not (types-equal-or-intersect (lvar-type x) (specifier-type 'single-float)))
+          (safe-ctype-for-single-coercion-p (lvar-type y)))
+      `(,(lvar-fun-name (basic-combination-fun node))
+         x (float y x))
+      (give-up-ir1-transform)))
 
 (dolist (x '(+ * / -))
   (%deftransform x '(function (rational float) *) #'float-contagion-arg1)
@@ -262,6 +334,68 @@
   (%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)
+    (multiple-value-bind (significand exponent sign)
+        ;; Signals an error for NaNs and infinities.
+        (handler-case (integer-decode-float x)
+          (error () (return-from maybe-exact-reciprocal nil)))
+      (let ((expected (/ sign significand (expt 2 exponent))))
+        (let ((reciprocal (/ 1 x)))
+          (multiple-value-bind (significand exponent sign) (integer-decode-float reciprocal)
+            (when (eql expected (* sign significand (expt 2 exponent)))
+              reciprocal)))))))
+
+;;; 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 subsctraction 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
@@ -566,7 +700,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
@@ -577,8 +711,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))
@@ -602,9 +736,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)
@@ -1108,27 +1242,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))
@@ -1137,39 +1289,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)))))))
 
@@ -1183,7 +1377,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
@@ -1204,8 +1398,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
@@ -1300,15 +1494,48 @@
   (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 (lvar-derived-type result)))
+                    (if (or (not y)
+                            (and (constant-lvar-p y) (= 1 (lvar-value y))))
+                        (if (values-type-p result-type)
+                            `(let ((res (,',unary x)))
+                               (values res (- x (,',coerce res))))
+                            `(let ((res (,',unary x)))
+                               ;; Dummy secondary value!
+                               (values res x)))
+                        (if (values-type-p result-type)
+                            `(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)))