What I had in mind is maybe less drastic then what you all are
purposing. I was only thinking about a function that would
evaluate a Maxima expression using native Lisp complex numbers
and return a float; it wouldn't change the representation of
any Maxima expression. Here is a demo of a first pass at
such a function
(C1) load("c:/cfloat.lisp");
Loading c:/cfloat.lisp
Finished loading c:/cfloat.lisp
(D1) c:/cfloat.lisp
(C2) z : [cos(1+%i), acos(1+%i), exp(2/%i), sqrt(3-%i)]$
(C3) cfloat(z);
(D3) [0.83373002513115 - 0.98889770576287 %I,
0.90455689430238 - 1.061275061905036 %I,
- 0.90929742682568 %I - 0.41614683654714,
1.755317301824428 - 0.28484878459314 %I]
(C4) xfloat(x) := scanmap('float, rectform(x))$
/* xfloat(z) gives the same value.*/
(C5) xfloat(z);
(D5) [0.83373002513115 - 0.98889770576287 %I,
0.90455689430238 - 1.061275061905036 %I,
- 0.90929742682568 %I - 0.41614683654714,
1.755317301824428 - 0.28484878459314 %I]
/* cfloat is usually faster than is xfloat. */
(C6) showtime : all$
Evaluation took 00.20 seconds (0.00 elapsed)
(C7) z : sum((1+%i/5)^k,k,0,199)$
Evaluation took 1.33 seconds (1.33 elapsed)
(C8) cfloat(z);
Evaluation took 0.25 seconds (0.25 elapsed)
(D8) 57.44462247189882 %I + 247.0188380491944
/* Let's check the value. */
(C9) cfloat(((1+%i/5)^200 - 1)/(%i/5));
Evaluation took 0.00 seconds (0.00 elapsed)
(D9) 57.44462247189757 %I + 247.0188380491941
/* Now evaluate with xfloat; we see that xfloat is much slower. In
this example, cfloat is more accurate (as verified via a bfloat
calculation); this won't always be true. */
(C10) xfloat(z);
Evaluation took 63.00 seconds (63.00 elapsed)
(D10) 57.44462502838694 %I + 247.0188405630827
/* cfloat is (supposed!) to give up when expression doesn't evaluate
to a float. */
(C11) cfloat(x + 5/6);
(sub)expression x doesn't evaluate to a float
-- an error. Quitting. To debug this try DEBUGMODE(TRUE);
Note: The above was done using maxima 5.5 under NT; it also worked
for maxima 5.9 under NT; however, under 5.9 xmaxima reports bogus
execution times.
gcl and cfloat have the share the same "weirdness":
(C15) cfloat((-1)^0.5);
Evaluation took 00.20 seconds (00.20 elapsed)
(D15) %I + 6.123031769111886E-17
(C16) to_lisp();
Type (run) to restart
MAXIMA>(sqrt -1)
#C(6.1230317691118863E-17 1.0)
MAXIMA>
Possible improvements and caveats to my code:
1. Look for %e ^ x and evaluate it as (exp x) instead of
(expt 2.718... x). Similarly, we could look for x ^ 0.5 and
use evaluate it as (sqrt x).
2. I'm sure I've missed a bunch of functions; the function
float-name probably doesn't transform everything it ought to.
3. For fun, we could (optionally) to use compensated addition.
-----------my minimally tested code ----------------------
;; Author: Barton Willis
;; University of Nebraska at Kearney (aka UNK, aka ...)
;; Warning: this code has only been minimally tested. Don't
;; trust it for anything important.
(defun $cfloat (e)
(cond ((or ($listp e) ($matrixp e))
(let ((op (car e)))
(cons op (mapcar #'$cfloat (cdr e)))))
(t
(let ((z (cfloat e nil)))
(cond ((complexp z)
(add (realpart z) (mul '$%i (imagpart z))))
(t
z))))))
(defun cfloat (e &optional fun)
(cond ((numberp e)
(float e 1.0L0))
((eq e '$%i)
#c(0.0L0 1.0L0))
((or (complexp e) (floatp e)) e)
((and (symbolp e) (mget e '$numer)))
(($bfloatp e)
(fp2flo e))
(($atom e)
(merror "(sub)expression ~:M doesn't evaluate to a float" e))
((eq (caar e) 'rat)
(fpcofrat e))
((setq fun (float-name (caar e)))
(setq e (mapcar #'cfloat (cdr e)))
(cond ((or (eq fun '+) (eq fun '*))
(reduce fun e))
(t
(apply fun e))))
(t
(merror "(sub)expression ~:M doesn't evaluate to a float" e))))
(declaim (inline float-name))
(defun float-name (f)
(cond ((eq f 'mplus)
'+)
((eq f 'mtimes)
'*)
((eq f 'mexpt)
'expt)
((eq (getchar f 1) '%)
(intern (string-left-trim "%" (symbol-name f)) 'lisp))
(t
nil)))
Regards,
Barton