Hello Richard,
thank you for your suggestions.
I do the work on Maxima just for fun. I have studied physics, but I am no longer
working as a physcist since a long time. In the earlier times of computers I had
done a lot of programming work in languages like Fortran, Pascal, C, C++.
Lisp is very new for me and I like it. But I do not know all the more
shophisticated approaches to use the language as optimal as possible. Most of
the language I have learned looking at the code of Maxima.
So the more directly bigfloat implementation is for me just a step more to get
experience in programming Lisp and Maxima.
Today, I have implemented a routine fpz-beta-incomplete. This routine does the
complex bigfloat evaluation more directly too. I have used the bigfloat code
presented by David Billinghurst in a post from January 15, 2006 "Code for
extended precision complex floating point arithmetic". I have only added one
routine. The code works wonderful. The gain in effiency is about 50% in
comparsion with an algorithm using only add, mul, ...
My question:
Should I add the complex bigfloat to the file float.lisp.
Next, I could commit an implementation of the Incomplete Beta function using the
more directly bigfloat implementations.
That is the code for complex bigfloat evaluation of Incomplete Beta:
*********************
(defun fpz-beta-incomplete (a b z)
(let* ((beta-maxit *beta-incomplete-maxit*)
(fpzero (cdr ($bfloat bigfloatzero)))
(fpone (cdr ($bfloat bigfloatone)))
(fptwo (fpplus fpone fpone))
(fpzzero (fpz fpzero fpzero))
(fpzone (fpz fpone fpzero))
(fpzmone (fpz (fpminus fpone) fpzero))
(fpztwo (fpz fptwo fpzero)))
(cond
((fpgreaterp
(fpz-re z)
(fpz-re
(fpz/
(fpz+ a fpzone)
(fpz+ a (fpz+ b fpztwo)))))
(when *debug*
(format t "~& in condition fpgreaterp~%"))
(fpz-
(fpz/
(fpz*
(bigfloattofpz (simplify (list '(%gamma) (fpztobigfloat a))))
(bigfloattofpz (simplify (list '(%gamma) (fpztobigfloat b)))))
(bigfloattofpz
(simplify (list '(%gamma) (fpztobigfloat (fpz+ a b))))))
(fpz-beta-incomplete b a (fpz- fpzone z))))
(t
(let* ((beta-eps (cdr ($bfloat (power ($bfloat 10.0) (- $fpprec)))))
(beta-min (fptimes* beta-eps beta-eps))
(fpz-beta-min (fpz beta-min fpzero))
(ab (fpz+ a b))
(ap (fpz+ a fpzone))
(am (fpz- a fpzone))
(c fpzone)
(d (fpz- fpzone (fpz/ (fpz* z ab) ap)))
(d (if (fplessp (fpzabs d) beta-min) fpz-beta-min d))
(d (fpz/ fpzone d))
(h d)
(aa fpzone)
(m2 fpzzero))
(do ((m 1 (+ m 1)))
((> m beta-maxit)
(merror "Continued fractions failed in `beta_incomplete'"))
(setq m2 (bigfloattofpz (+ m m)))
(setq aa (fpz/
(fpz*
(bigfloattofpz m)
(fpz* z (fpz- b (bigfloattofpz m))))
(fpz* (fpz+ am m2) (fpz+ a m2))))
(setq d (fpz+ fpzone (fpz* aa d)))
(when (fplessp (fpzabs d) beta-min) (setq d fpz-beta-min))
(setq c (fpz+ fpzone (fpz/ aa c)))
(when (fplessp (fpzabs c) beta-min) (setq c fpz-beta-min))
(setq d (fpz/ fpzone d))
(setq h (fpz* h (fpz* d c)))
(setq aa (fpz/
(fpz*
fpzmone
(fpz*
(fpz+ a (bigfloattofpz m))
(fpz* (fpz+ ab (bigfloattofpz m)) z)))
(fpz*
(fpz+ a m2) (fpz+ ap m2))))
(setq d (fpz+ fpzone (fpz* aa d)))
(when (fplessp (fpzabs d) beta-min) (setq d fpz-beta-min))
(setq c (fpz+ fpzone (fpz/ aa c)))
(when (fplessp (fpzabs c) beta-min) (setq c fpz-beta-min))
(setq d (fpz/ fpzone d))
(setq de (fpz* d c))
(setq h (fpz* h de))
(when (fplessp (fpzabs (fpz- de fpzone)) beta-eps)
(return
(fpz/
(fpz*
h
(fpz*
(bigfloattofpz
(power (fpztobigfloat z) (fpztobigfloat a)))
(bigfloattofpz
(power
(fpztobigfloat (fpz- fpzone z)) (fpztobigfloat b)))))
a)))))))))
*********************
Dieter Kaiser