diff --git a/compiler/float-tran.lisp b/compiler/float-tran.lisp index 9c066ab3fa058e842ec1aa97a9bf674b29c9d893..c83ad49f7f45bc670b1219d26816ec22063e4feb 100644 --- a/compiler/float-tran.lisp +++ b/compiler/float-tran.lisp @@ -5,7 +5,7 @@ ;;; Carnegie Mellon University, and has been placed in the public domain. ;;; (ext:file-comment - "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/compiler/float-tran.lisp,v 1.121.2.1.2.1 2008/10/09 16:25:47 rtoy Exp $") + "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/compiler/float-tran.lisp,v 1.121.2.1.2.2 2008/10/10 04:08:40 rtoy Exp $") ;;; ;;; ********************************************************************** ;;; @@ -1452,7 +1452,8 @@ (ry (realpart y)) (iy (imagpart y))) (complex (- (* rx ry) (* ix iy)) - (+ (* rx iy) (* ix ry))))) + (+ (* rx iy) (* ix ry))))) + #+nil (deftransform / ((x y) ((complex ,type) (complex ,type)) * :policy (> speed space)) '(let* ((rx (realpart x)) @@ -1474,6 +1475,7 @@ (deftransform * ((z w) (real (complex ,type)) *) '(complex (* (realpart w) z) (* (imagpart w) z))) ;; Divide a complex by a real + #+nil (deftransform / ((w z) ((complex ,type) real) *) '(complex (/ (realpart w) z) (/ (imagpart w) z))) ;; Divide a real by a complex @@ -1526,6 +1528,64 @@ (complex (- (* rx ry) (* ix iy)) (+ (* rx iy) (* ix ry))))) +(deftransform / ((x y) ((complex single-float) (complex single-float)) * + :policy (> speed space)) + '(let* ((rx (realpart x)) + (ix (imagpart x)) + (ry (realpart y)) + (iy (imagpart y))) + (if (> (abs ry) (abs iy)) + (let* ((r (/ iy ry)) + (dn (+ ry (* r iy)))) + (complex (/ (+ rx (* ix r)) dn) + (/ (- ix (* rx r)) dn))) + (let* ((r (/ ry iy)) + (dn (+ iy (* r ry)))) + (complex (/ (+ (* rx r) ix) dn) + (/ (- (* ix r) rx) dn)))))) + +(deftransform / ((x y) ((complex double-float) (complex double-float)) * + :policy (> speed space)) + ;; Divide a complex by a complex + + ;; Here's how we do a complex division + ;; + ;; Compute (xr + i*xi)/(yr + i*yi) + ;; + ;; Assume |yi| < |yr|. Then + ;; + ;; (xr + i*xi) (xr + i*xi) + ;; ----------- = ----------------- + ;; (yr + i*yi) yr*(1 + i*(yi/yr)) + ;; + ;; (xr + i*xi)*(1 - i*(yi/yr)) + ;; = --------------------------- + ;; yr*(1 + (yi/yr)^2) + ;; We do the similar thing when |yi| > |yr|. The result is + ;; + ;; + ;; (xr + i*xi) (xr + i*xi) + ;; ----------- = ----------------- + ;; (yr + i*yi) yi*((yr/yi) + i) + ;; + ;; (xr + i*xi)*((yr/yi) - i) + ;; = ------------------------- + ;; yi*((yr/yi)^2 + 1) + '(let* ((ry (realpart y)) + (iy (imagpart y))) + (if (> (abs ry) (abs iy)) + (let* ((r (/ iy ry)) + (dn (+ ry (* r iy)))) + (/ (* x (complex 1d0 r)) + dn)) + (let* ((r (/ ry iy)) + (dn (+ iy (* r ry)))) + (/ (* x (complex r -1d0)) + dn))))) + +(deftransform / ((w z) ((complex single-float) real) *) + '(complex (/ (realpart w) z) (/ (imagpart w) z))) + #+(and (not double-double) complex-fp-vops) (macrolet ((frob (type) diff --git a/compiler/x86/float-sse2.lisp b/compiler/x86/float-sse2.lisp index 539030e7e316a62243afb185c8a2042e327ecadd..bd32b0ebd1aea67a7b29870b6b428dcf81b80ef6 100644 --- a/compiler/x86/float-sse2.lisp +++ b/compiler/x86/float-sse2.lisp @@ -7,7 +7,7 @@ ;;; Scott Fahlman or slisp-group@cs.cmu.edu. ;;; (ext:file-comment - "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/compiler/x86/float-sse2.lisp,v 1.1.2.8.2.5 2008/10/10 03:11:33 rtoy Exp $") + "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/compiler/x86/float-sse2.lisp,v 1.1.2.8.2.6 2008/10/10 04:08:40 rtoy Exp $") ;;; ;;; ********************************************************************** ;;; @@ -2047,3 +2047,19 @@ (inst shufpd t2 t2 1) ; t2 = a*d|-b*d (inst addpd t2 t1) ; t2 = a*d+b*c | a*c-b*d (inst movapd r t2))) + +;; Divide a complex by a double-float +(define-vop (//complex-double-float/double-float) + (:translate /) + (:args (x :scs (complex-double-reg)) + (y :scs (double-reg))) + (:arg-types complex-double-float double-float) + (:results (r :scs (complex-double-reg))) + (:result-types complex-double-float) + (:policy :fast-safe) + (:temporary (:scs (complex-double-reg)) t0) + (:generator 1 + (inst movq t0 y) ; t0 = 0|y + (inst unpcklpd t0 t0) ; t0 = y|y + (inst movapd r x) ; r = b|a (x = a+b*i) + (inst divpd r t0))) ; r = b/y|a/y