;;; -*- Mode:LISP; Package:NEW-MATH; Base:10; Readtable:CL -*-


;-----------------------------------------------------------------------------
;   Realpart of complex
;-----------------------------------------------------------------------------
(defun realpart (x)
  (prims:dispatch vinc::%%data-type x
		  (($$dtp-fixnum $$dtp-bignum $$dtp-short-float $$dtp-single-float
				 $$dtp-double-float $$dtp-rational) x)
		  ($$dtp-complex (progn
				   (hw:vma-start-read-vma-boxed-md-boxed x)
				   (hw:read-md)))
     (t (trap::illop "Can't take realpart of that"))))

;-----------------------------------------------------------------------------
;   Imagpart of complex
;-----------------------------------------------------------------------------
(defun imagpart (x)
  (prims:dispatch vinc::%%data-type x
		  (($$dtp-fixnum $$dtp-bignum $$dtp-short-float $$dtp-single-float
				 $$dtp-double-float $$dtp-rational) 0)
		  ($$dtp-complex (progn
				   (hw:vma-start-read-cdr-vma-boxed-md-boxed x)
				   (hw:read-md)))
     (t (trap::illop "Can't take imagpart of that"))))

;-----------------------------------------------------------------------------
;   make canonical complex
;-----------------------------------------------------------------------------
(defun complex (real imag)
  (when (or (not (numberp real)) (complexp real)
	    (not (numberp imag)) (complexp imag))
    (li:error "Bad arguments to COMPLEX"))
  (multiple-value-bind (real imag)
      (generic-math-type-coercer real imag)
    (if (and (not (floatp imag)) (zerop imag))
	(values real (multiple-value-bind (num status)
			 (test-generic real)
		       status))
      (let ((result (hw:dpb-boxed vinc:$$dtp-complex vinc::%%data-type (cons:cons real imag))))
	(multiple-value-bind (ignore status)
	    (test-complex result)
	  (values result status))))))

;-----------------------------------------------------------------------------
;   add complex
;-----------------------------------------------------------------------------
(defun add-complex (x y)
  (complex
    (add-generic (realpart x) (realpart y))
    (add-generic (imagpart x) (imagpart y))))
;-----------------------------------------------------------------------------
;   subtract complex
;-----------------------------------------------------------------------------
(defun subtract-complex (x y)
  (complex
    (subtract-generic (realpart x) (realpart y))
    (subtract-generic (imagpart x) (imagpart y))))
;-----------------------------------------------------------------------------
;   multiply complex
;-----------------------------------------------------------------------------
(defun multiply-complex (x y)
  (let* ((x-real (realpart x))
	 (y-real (realpart y))
	 (x-imag (imagpart x))
	 (y-imag (imagpart y)))
    (complex
      (subtract-generic
	(multiply-generic x-real y-real)
	(multiply-generic x-imag y-imag))
      (add-generic
	(multiply-generic x-real y-imag)
	(multiply-generic x-imag y-real)))))
;-----------------------------------------------------------------------------
;   divide complex
;-----------------------------------------------------------------------------
(defun divide-complex (x y)
  (let* ((x-real (realpart x))
	 (y-real (realpart y))
	 (x-imag (imagpart x))
	 (y-imag (imagpart y))
	 (divisor (add-generic
		    (multiply-generic y-real y-real)
		    (multiply-generic y-imag y-imag))))
    (complex
      (divide-generic
	(add-generic
	  (multiply-generic x-real y-real)
	  (multiply-generic x-imag y-imag))
	divisor)
      (divide-generic
	(subtract-generic
	  (multiply-generic x-imag y-real)
	  (multiply-generic x-real y-imag))
	divisor))))

;-----------------------------------------------------------------------------
;   negate complex
;-----------------------------------------------------------------------------
(defun negate-complex (x)
  (complex
    (- (realpart x))
    (- (imagpart x))))

;-----------------------------------------------------------------------------
;   compare complex
;-----------------------------------------------------------------------------
(defun compare-complex (x y) ;only equality is tested for!
  (if (and (= (realpart x) (realpart y))
	   (= (imagpart x) (imagpart y)))
      (values 0 #x040000)
    (values 1 0)))

;-----------------------------------------------------------------------------
;   test complex
;-----------------------------------------------------------------------------
(defun test-complex (x) ;only test for zerop
  (if (and (zerop (realpart x)) (zerop (imagpart x)))
      (values 0 #x040000)
    (values 1 0)))

;-----------------------------------------------------------------------------
;   complex conjugate
;-----------------------------------------------------------------------------
(defun conjugate (number)
  (complex
    (realpart number)
    (- (imagpart number))))