Skip to content
Snippets Groups Projects
Commit bf294fb3 authored by rtoy's avatar rtoy
Browse files

x86/float-sse2.lisp:

o Add vop to divide a complex by a real using packed operations.

float-tran.lisp:
o New transform for /.  Take advantage of packed complex multiply and
  packed division of a complex by a real.
o Remove deftransform for dividing a complex double-float by a
  double-float.

These give a moderate speedup for complex double-float divison.  Still
seems much slower than x87.
parent c8d26008
No related branches found
No related tags found
No related merge requests found
......@@ -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)
......
......@@ -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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment