#|
Original
(COMPLEX IDENTITY (COMPLEX)) (single/double) 905 / 906
(COMPLEX - (COMPLEX)) (single/double) 2315 / 2304
(COMPLEX + (COMPLEX REAL)) (single/double) 1926 / 1921
(COMPLEX + (REAL COMPLEX)) (single/double) 1925 / 1920
(COMPLEX + (COMPLEX COMPLEX)) (single/double) 2434 / 2434
(COMPLEX - (COMPLEX REAL)) (single/double) 1926 / 1921
(COMPLEX - (REAL COMPLEX)) (single/double) 2441 / 2436
(COMPLEX - (COMPLEX COMPLEX)) (single/double) 2436 / 2436
(COMPLEX * (COMPLEX REAL)) (single/double) 1990 / 1928
(COMPLEX * (REAL COMPLEX)) (single/double) 1925 / 2027
(COMPLEX * (COMPLEX COMPLEX)) (single/double) 2952 / 2952
(COMPLEX / (COMPLEX (REAL 1))) (single/double) 7421 / 9842
(COMPLEX / (REAL (COMPLEX #C(1 1)))) (single/double) 14351 / 17680
(COMPLEX / ((COMPLEX #C(1 1)) (COMPLEX #C(-2 -2)))) (single/double) 12810 / 15565
(COMPLEX CONJUGATE (COMPLEX)) (single/double) 1422 / 1677
(T = (COMPLEX REAL)) (single/double) 2189 / 1932
(T = (REAL COMPLEX)) (single/double) 2189 / 1951
(T = (COMPLEX COMPLEX)) (single/double) 2189 / 2189

With VOPs:
(COMPLEX IDENTITY (COMPLEX)) (single/double) 653 / 771
(COMPLEX - (COMPLEX)) (single/double) 1421 / 1205
(COMPLEX + (COMPLEX REAL)) (single/double) 1022 / 915
(COMPLEX + (REAL COMPLEX)) (single/double) 1025 / 915
(COMPLEX + (COMPLEX COMPLEX)) (single/double) 910 / 923
(COMPLEX - (COMPLEX REAL)) (single/double) 1022 / 915
(COMPLEX - (REAL COMPLEX)) (single/double) 1025 / 934
(COMPLEX - (COMPLEX COMPLEX)) (single/double) 934 / 928
(COMPLEX * (COMPLEX REAL)) (single/double) 1143 / 1005
(COMPLEX * (REAL COMPLEX)) (single/double) 1264 / 1187
(COMPLEX * (COMPLEX COMPLEX)) (single/double) 2023 / 2190
(COMPLEX / (COMPLEX (REAL 1))) (single/double) 4857 / 4599
(COMPLEX / (REAL (COMPLEX #C(1 1)))) (single/double) 14871 / 17511
(COMPLEX / ((COMPLEX #C(1 1)) (COMPLEX #C(-2 -2)))) (single/double) 12815 / 14094
(COMPLEX CONJUGATE (COMPLEX)) (single/double) 1086 / 1275
(T = (COMPLEX REAL)) (single/double) 1667 / 1425 ; EQL is even better
(T = (REAL COMPLEX)) (single/double) 1585 / 1426
(T = (COMPLEX COMPLEX)) (single/double) 1473 / 1473

Few VOPs:
(COMPLEX IDENTITY (COMPLEX)) (single/double) 653 / 771
(COMPLEX - (COMPLEX)) (single/double) 1421 / 1205
(COMPLEX + (COMPLEX REAL)) (single/double) 2314 / 1509
(COMPLEX + (REAL COMPLEX)) (single/double) 3201 / 2430
(COMPLEX + (COMPLEX COMPLEX)) (single/double) 935 / 910
(COMPLEX - (COMPLEX REAL)) (single/double) 2313 / 1508
(COMPLEX - (REAL COMPLEX)) (single/double) 3201 / 2430
(COMPLEX - (COMPLEX COMPLEX)) (single/double) 935 / 910
(COMPLEX * (COMPLEX REAL)) (single/double) 1415 / 1416
(COMPLEX * (REAL COMPLEX)) (single/double) 2423 / 2424
(COMPLEX * (COMPLEX COMPLEX)) (single/double) 3461 / 3460
(COMPLEX / (COMPLEX (REAL 1))) (single/double) 4857 / 4597
(COMPLEX / (REAL (COMPLEX #C(1 1)))) (single/double) 12814 / 12706
(COMPLEX / ((COMPLEX #C(1 1)) (COMPLEX #C(-2 -2)))) (single/double) 15239 / 16009
(COMPLEX CONJUGATE (COMPLEX)) (single/double) 1086 / 1275
(T = (COMPLEX REAL)) (single/double) 1667 / 1426
(T = (REAL COMPLEX)) (single/double) 1586 / 1426
(T = (COMPLEX COMPLEX)) (single/double) 1461 / 1474
|#

;; for some reason the sequence below isn't open coded...

(in-package "SB-C")
(deftransform / ((x y) (single-float (complex single-float)) *)
  `(/ (complex x ,(coerce 0 'single-float)) y))
(deftransform / ((x y) (double-float (complex double-float)) *)
  `(/ (complex x ,(coerce 0 'double-float)) y))
;;; Define some transforms for complex operations. We do this in lieu
;;; of complex operation VOPs.
(macrolet ((frob (type)
             `(progn
                #-sb-xc-host
                (deftransform complex ((r) (,type))
                  `(complex r ,(coerce 0 ',type)))
               ;; negation
                #-x86-64
               (deftransform %negate ((z) ((complex ,type)) *)
                 '(complex (%negate (realpart z)) (%negate (imagpart z))))
               ;; complex addition and subtraction
               #-x86-64
               (deftransform + ((w z) ((complex ,type) (complex ,type)) *)
                 '(complex (+ (realpart w) (realpart z))
                           (+ (imagpart w) (imagpart z))))
               #-x86-64
               (deftransform - ((w z) ((complex ,type) (complex ,type)) *)
                 '(complex (- (realpart w) (realpart z))
                           (- (imagpart w) (imagpart z))))
               ;; Add and subtract a complex and a real.
               (deftransform + ((w z) ((complex ,type) real) *)
                 #-x86-64
                 '(complex (+ (realpart w) z) (imagpart w))
                 #+x86-64
                 `(+ w (complex z)))
               (deftransform + ((z w) (real (complex ,type)) *)
                 #-x86-64
                 '(complex (+ (realpart w) z) (imagpart w))
                 #+x86-64
                 `(+ (complex z) w))
               ;; Add and subtract a real and a complex number.
               (deftransform - ((w z) ((complex ,type) real) *)
                 #-x86-64
                 '(complex (- (realpart w) z) (imagpart w))
                 #+x86-64
                 `(- w (complex z)))
               (deftransform - ((z w) (real (complex ,type)) *)
                 #-x86-64
                 '(complex (- z (realpart w)) (- (imagpart w)))
                 #+x86-64
                 `(- (complex z) w))
               ;; Multiply and divide two complex numbers.
               (deftransform * ((x y) ((complex ,type) (complex ,type)) *)
                 '(let* ((rx (realpart x))
                         (ix (imagpart x))
                         (ry (realpart y))
                         (iy (imagpart y)))
                   (declare (ignorable ry iy))
                   #+x86-64
                   (+ (%pairwise-* (complex rx rx) y)
                      (%pairwise-* (complex ix ix) (%negate-real (%complex-swap y))))
                   #-x86-64
                    (complex (- (* rx ry) (* ix iy))
                             (+ (* rx iy) (* ix ry)))))
               (deftransform / ((x y) (,type (complex ,type)) *)
                 `(/ (complex x) y))
               (deftransform / ((x y) ((complex ,type) (complex ,type)) *)
                 '(let* ((rx (realpart x))
                         (ix (imagpart x))
                         (ry (realpart y))
                         (iy (imagpart y))
                         (swap (%complex-swap (%negate-real x))))
                   (declare (ignorable rx ix))
                   (multiple-value-call #'/
                     (if (> (abs ry) (abs iy))
                         (let* ((r (/ iy ry))
                                (dn (* ry (+ 1 (* r r)))))
                           (values (+ x (* swap r)) dn))
                         (let* ((r (/ ry iy))
                                (dn (* iy (+ 1 (* r r)))))
                           (values (+ (* x r) swap) dn))))))
               ;; Multiply a complex by a real or vice versa.
               (deftransform * ((w z) ((complex ,type) real) *)
                 #-x86-64
                 '(complex (* (realpart w) z) (* (imagpart w) z))
                 #+x86-64
                 '(%pairwise-* w (complex z z)))
               (deftransform * ((z w) (real (complex ,type)) *)
                 #-x86-64
                 '(complex (* (realpart w) z) (* (imagpart w) z))
                 #+x86-64
                 '(%pairwise-* (complex z z) w))
               ;; Divide a complex by a real.
               #-x86-64
               (deftransform / ((w z) ((complex ,type) real) *)
                 '(complex (/ (realpart w) z) (/ (imagpart w) z)))
               ;; conjugate of complex number
               #-x86-64
               (deftransform conjugate ((z) ((complex ,type)) *)
                 '(complex (realpart z) (- (imagpart z))))
               ;; CIS
               (deftransform cis ((z) ((,type)) *)
                 '(complex (cos z) (sin z)))
               ;; comparison
               #-x86-64
               (deftransform = ((w z) ((complex ,type) (complex ,type)) *)
                 '(and (= (realpart w) (realpart z))
                       (= (imagpart w) (imagpart z))))
               #-x86-64
               (deftransform = ((w z) ((complex ,type) real) *)
                 '(and (= (realpart w) z) (zerop (imagpart w))))
               #-x86-64
               (deftransform = ((w z) (real (complex ,type)) *)
                 '(and (= (realpart z) w) (zerop (imagpart z)))))))

  (frob single-float)
  (frob double-float))
(in-package "CL-USER")

(defun bench (n fun)
  (declare (type unsigned-byte n))
  (unless (functionp fun)
    (setf fun (fdefinition fun)))
  (aref (sort (map-into (make-array n :element-type `(unsigned-byte ,sb-vm:n-word-bits))
                        (lambda ()
                          (declare (type function fun))
                          (nth-value
                           1
                           (sb-vm::with-cycle-counter
                             (funcall fun)))))
              #'<)
        (truncate n 2)))

(defmacro test-fun (fun ret (&rest args) &key (n 100000) (m 256) (float-type 'double-float))
  (let* ((inits '())
         (types (prog1
                    (mapcar (lambda (type)
                              (push (and (consp type)
                                         (prog1 (second type)
                                           (setf type (car type))))
                                    inits)
                              (case type
                                (real
                                   float-type)
                                (complex
                                   `(complex ,float-type))
                                (t type)))
                            (cons ret args))
                  (setf inits (nreverse inits))))
         (vars (loop for x in types
                     collect (gensym "VEC"))))
    `(let ,(mapcar (lambda (var type init)
                     `(,var (make-array ,m :element-type ',type
                                        ,@(and init
                                               `(:initial-element ,(coerce init type))))))
            vars types inits)
       (declare ,@(mapcar (lambda (var type)
                            `(type (simple-array ,type (,m)) ,var))
                          vars types))
       (bench ,n (lambda ()
                   (declare (optimize speed))
                   (map-into ,(first vars) #',fun ,@(rest vars)))))))

(defun run ()
  (macrolet
      ((run (&body settings)
         `(progn
            ,@(loop for (ret fun . args) in settings
                    collect
                    `(format t "~S (single/double) ~A / ~A~%"
                             '(,ret ,fun ,args)
                              (test-fun ,fun ,ret ,args
                                   :float-type single-float)
                              (test-fun ,fun ,ret ,args
                                   :float-type double-float))))))
    (run (complex identity complex)
         (complex - complex)
         (complex + complex real)
         (complex + real complex)
         (complex + complex complex)
         (complex - complex real)
         (complex - real complex)
         (complex - complex complex)
         (complex * complex real)
         (complex * real complex)
         (complex * complex complex)
         (complex / complex (real 1))
         (complex / real (complex #c(1 1)))
         (complex / (complex #c(1 1)) (complex #c(-2 -2)))
         (complex conjugate complex)
         (t       = complex real)
         (t       = real complex)
         (t       = complex complex))))