[RFC] Code for extended precision complex floating point arithmetic

from the perspective of the author of the macsyma bigfloat stuff, this
looks ok, relative to programming circa 1974.  It is possible that one
should write this differently given Common Lisp and CLOS.
 I would organize it differently myself.  See also comments about MPFR too.

For example, instead of using lists for <fraction, exponent> pairs,
use structures.  The names here are too long but ..  something like
(defstruct complex-bigfloat  ( real-part bfzero :type bigfloat) ( imag-part 
bfzero  :type bigfloat)  precision)

Instead of defining complex bigfloat +  as fpz+,   use what amounts to
(defmethod two-arg-+ ( (x complex-bigfloat)(y complex-bigfloat))    ....
which also allows you to defineb
(defmethod two-arg-+ ((x complex-bigfloat) (y  bigfloat)) ...  and
(defmethod two-arg-+ ((x bigfloat) (y complex-bigfloat)) ...    and
(defmethod two-arg-+ ((x complex-bigfloat)( y number))  ...  for all common 
lisp number types.

Oh, you define generic  +   in terms of two-arg-+,  so you can also do (+ a 
b c)   for complex bigfloats a b c.

What does this get you:
  More versatility
 Probably slower speed

The same program that computes 2F1 for complex bigfloat, may also work for
double-float [unchanged except for different epsilon parameters etc]. Though 
if
you know you only need double-float there are probably shortcuts.

 The way to use more bits for division... bump up the floating precision
if you really need it.. The division should be correctly rounded.

Also, if you need names for these operations, look at the GMP names or 
perhaps
the float package on top of GMP  http://www.mpfr.org/

This has probably made your life more complicated, if you care to follow it 
all.

but if you do, I can give you more details and code for generic arithmetic, 
and
it could possibly be a model for other improvements to Maxima.

Also, what if we just said we would do a really bang-up job of double-double
(i.e. quad precision), and imported a really good library for this.  Would 
we
still have much call for bigfloats?
Or could we kind of ooze over from double to quad to slow bigfloats?

Or should we really be importing the MPFR library?

I found this, which is claimed to be generic code for hypergeometric
series...

http://cvs.rpm.org/viewcvs/genius/mpfr/generic.c?rev=1.2


Sorry for complicating things..
RJF



----- Original Message ----- 
From: CALCRTS" 
<david.billinghurst at comalco.riotinto.com.au>
To: <maxima at math.utexas.edu>
Sent: Sunday, January 15, 2006 3:53 PM
Subject: [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.
>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>