;; complexfloat.lisp
;;
;; Routines for extended precision complex floating point numbers
;; See float.lisp

;; Complex FP number x = xr + %i*xi is represented by (list xr xi)

(in-package :maxima)

;; define the data structure for a complex FP number
;; Complex FP number x = xr + %i*xi is represented by (list xr xi)
(defun fpz-re (z) (car z))
(defun fpz-im (z) (cadr z))
(defun fpz (x y) `(,x ,y))

;; z = x + y
(defun fpz+ (x y)
  (fpz (fpplus (fpz-re x) (fpz-re y)) 
       (fpplus (fpz-im x) (fpz-im y))))

;; z = x - y
(defun fpz- (x y)
  (fpz (fpdifference (fpz-re x) (fpz-re y))
       (fpdifference (fpz-im x) (fpz-im y))))

;; z = x * y   FIXME: Probably want a guard digit or two
(defun fpz* (x y)
  (fpz (fpdifference (fptimes* (fpz-re x) (fpz-re y)) 
		     (fptimes* (fpz-im x) (fpz-im y)))
       (fpplus (fptimes* (fpz-re x) (fpz-im y))
	       (fptimes* (fpz-im x) (fpz-re y)))))

;; z = x / y  FIXME: Probably want a guard digit or two
(defun fpz/ (x y)
  (cond ((and (equal (fpz-re y) '(0 0)) (equal (fpz-im y) '(0 0)))
	 (merror "complex division by zero"))
	((equal (fpz-im y) '(0 0))
	 (fpz (fpquotient (fpz-re x) (fpz-re y)) 
	      (fpquotient (fpz-im x) (fpz-re y))))
	((equal (fpz-re y) '(0 0))
	 (fpz (fpquotient (fpz-im x) (fpz-im y)) 
	       (fpquotient (fpminus (fpz-re x)) (fpz-im y))))
	(t (let
	     ((denominator (fpplus (fptimes* (fpz-re y) (fpz-re y))
				   (fptimes* (fpz-im y) (fpz-im y))))
	      (re-numerator (fpplus (fptimes* (fpz-re x) (fpz-re y))
				    (fptimes* (fpz-im x) (fpz-im y))))
	      (im-numerator (fpdifference (fptimes* (fpz-im x) (fpz-re y))
					  (fptimes* (fpz-re x) (fpz-im y)))))
	     (fpz (fpquotient re-numerator denominator)
		  (fpquotient im-numerator denominator))))))

;; z = abs(x) for complex FP x 
;; Note: fproot takes a bigfloat
(defun fpzabs (x)
  (fproot (bcons (fpplus (fptimes* (fpz-re x) (fpz-re x)) 
			 (fptimes* (fpz-im x) (fpz-im x)))) 
	  2))


;; Some utilities for testing

;; Convert a maxima number to complex FP
(defun floattofpz (z)
  (let (($float2bf t))
    (fpz (intofp ($realpart z)) (intofp ($imagpart z)))))

;; Convert FP to float
(defun fptofloat (x) (fp2flo (bcons x)))

;; Convert FP complex to maxima float
;; Return real or im
(defun fpztofloat (z)
  (cond ((equal (fpz-im z) '(0 0))
	   (fptofloat (fpz-re z)))
	 ((equal (fpz-re z) '(0 0))
	   `((mtimes simp) ,(fptofloat (fpz-im z)) $%i))
         (t `((mplus simp) ,(fptofloat (fpz-re z)) 
	      ((mtimes simp) ,(fptofloat (fpz-im z)) $%i)))))

;; Convert FP complex to maxima complex bigfloat
(defun fpztobigfloat (z)
  (cond ((equal (fpz-im z) '(0 0))
	   (bcons (fpz-re z)))
	 ((equal (fpz-re z) '(0 0))
	   `((mtimes simp) ,(bcons (fpz-im z)) $%i))
         (t `((mplus simp) ,(bcons (fpz-re z)) 
	      ((mtimes simp) ,(bcons (fpz-im z)) $%i)))))