[RFC] Code for extended precision complex floating point arithmetic

I have been writing some code to evaluate the hypergeometric 
functions 1F1(a;b,z) and 2F1(a,b;c;z) for complex arguments.
I now have working routines 

I needed to evaluate compex power series in extended precision,
so I wrote the code below.  It works, but I'd like to 
discuss the structure.  (My Fortran heritage shows through). 

The code is layered on top of the routines in float.lisp.

A maxima bigfloat is 
;; Representation of a Bigfloat:  ((BIGFLOAT SIMP precision) mantissa exponent)
;; precision -- number of bits of precision in the mantissa.  
;; 		precision = (haulong mantissa)
;; mantissa -- a signed integer representing a fractional portion computed by
;; 	       fraction = (// mantissa (^ 2 precision)).
;; exponent -- a signed integer representing the scale of the number.
;; 	       The actual number represented is (f* fraction (^ 2 exponent)).

Many of the lower level in float.lisp routines work on "floating point numbers"
which are just (mantissa exponent) pairs using the current value of fpprec.
I have chosen to represent "complex floating point numbers" as 
 (list (mantissa exponent) (mantissa exponent))

We then have
  fpz     Construct complex floating point number
  fpz-re  Real part of complex floating point number
  fpz-im  Imaginary part of complex floating point number
  fpz+    + for complex floating point number
  fpz-    - for complex floating point number
  fpz*    * for complex floating point number
  fpz/    / for complex floating point number
  fpzabs  Absolute value 

I would really welcome some feedback, particularly on
  - the data structure
  - function names.  
  - the FIXME comments on the code regarding algorithms.  I need to 
    reread Kahan's paper "Branch Cuts for Complex Elementary Functions"
    for a start).

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

(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 (x y) "FP complex from two FP numbers" `(,x ,y))
(defun fpz-re (z) "Real part of FP number" (car z))
(defun fpz-im (z) "Imaginary part of FP number" (cadr z))

(defun fpz+ (&rest args)
  "Add FP complex numbers"
  (let ((l (length args)))
    (cond ((equal l 0) (floattofpz 0))
	  ((equal l 1) (car args))
          ((equal l 2) 
	   (let ((x (car args)) (y (cadr args)))
	     (fpz (fpplus (fpz-re x) (fpz-re y)) 
		  (fpplus (fpz-im x) (fpz-im y)))))
	  (t (reduce #'fpz+ args)))))

(defun fpz- (x &rest args)
  "Difference of two FP complex numbers"
  (let ((l (length args)))
    (cond ((equal l 0) 
	    (fpz (fpminus (fpz-re x)) (fpminus (fpz-im x))))
          ((equal l 1) 
	   (let ((y (car args)))
	     (fpz (fpdifference (fpz-re x) (fpz-re y)) 
		  (fpdifference (fpz-im x) (fpz-im y)))))
	  (t (fpz- x (reduce #'fpz+ args))))))

(defun fpz* (&rest args)
  "Multilpy FP complex numbers"
  (let ((l (length args)))
    (cond ((equal l 0) (floattofpz 1))
	  ((equal l 1) (car args))
          ((equal l 2) 
	   (let ((x (car args)) (y (cadr args)))
	     (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))))))
	  (t (reduce #'fpz* args)))))

;; z = x / y  FIXME: Are guard digits required?
(defun fpz/ (x &rest args)
  "Quotient of FP complex numbers"
  (let ((l (length args)))
    (cond ((equal l 0) 
	    (fpz/ (floattofpz 1) x))
	  ((equal l 1) 
	   (let ((y (car args)))
	   (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)))))))
	  (t (fpz/ x (reduce #'fpz* args))))))

;; z = abs(x)  FIXME: Is a more sophisticated algorithm needed?
;; Note: fproot takes a bigfloat
(defun fpzabs (x)
  "Absolute value of a complex FP number"
  (fproot (bcons (fpplus (fptimes* (fpz-re x) (fpz-re x)) 
			 (fptimes* (fpz-im x) (fpz-im x)))) 
	  2))		       

;; 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)))))

;; Exponent of an FP number.  Return large -ve number if zero
(defun fp-exponent (x) 
  (if (equal (car x) 0) most-negative-fixnum (cadr x)))

;; Given a complex FP number, return larger exponent
(defun fpz-exponent (z)
  (max (fp-exponent (fpz-re z)) (fp-exponent (fpz-im z))))


NOTICE
This e-mail and any attachments are private and confidential and may contain privileged information. If you are not an authorised recipient, the copying or distribution of this e-mail and any attachments is prohibited and you must not read, print or act in reliance on this e-mail or attachments.
This notice should not be removed.