I have finished some work on the Incomplete Beta function. This time I have
implemented the bigfloat code very directly.
With fpprec:4096 the algorithm below needs on my system 55 seconds for an answer
(without display). The more native algorithm using add, mul ,... needs 82
seconds.
The gain in efficiency is about 30%. Not more.
Is there something which could be done much better (faster)?
Dieter Kaiser
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Bigfloat algorithm for numerical evaluation of Incomplete Beta
;;; a, b, z are pairs (mantissa exponent)
;;; The result is of the form (mantissa exponent) too.
(defun fp-beta-incomplete (a b z)
(let* ((beta-maxit *beta-incomplete-maxit*)
(fpzero (cdr ($bfloat bigfloatzero)))
(fpone (cdr ($bfloat bigfloatone)))
(fptwo (fpplus fpone fpone)))
(cond
((fpgreaterp
z
(fpquotient
(fpplus a fpone)
(fpplus a (fpplus b fptwo))))
(fpdifference
(fpquotient
(fptimes*
(cdr (simplify (list '(%gamma) (bcons a))))
(cdr (simplify (list '(%gamma) (bcons b)))))
(cdr (simplify (list '(%gamma) (bcons (fpplus a b))))))
(fp-beta-incomplete b a (fpdifference fpone z))))
(t
(let* ((beta-eps (cdr ($bfloat (power ($bfloat 10.0) (- $fpprec)))))
(beta-min (fptimes* beta-eps beta-eps))
(ab (fpplus a b))
(ap (fpplus a fpone))
(am (fpdifference a fpone))
(c fpone)
(d (fpdifference fpone (fpquotient (fptimes* z ab) ap)))
(d (if (fplessp d beta-min) beta-min d))
(d (fpquotient fpone d))
(h d)
(aa fpone)
(m2 fpzero))
(do* ((m 1 (+ m 1)))
((> m beta-maxit)
(merror "Continued fractions failed in `beta_incomplete'"))
(setq m2 (intofp (+ m m)))
(setq aa (fpquotient
(fptimes*
(intofp m)
(fptimes* z (fpdifference b (intofp m))))
(fptimes* (fpplus am m2) (fpplus a m2))))
(setq d (fpplus fpone (fptimes* aa d)))
(when (fplessp (fpabs d) beta-min) (setq d beta-min))
(setq c (fpplus fpone (fpquotient aa c)))
(when (fplessp (fpabs c) beta-min) (setq c beta-min))
(setq d (fpquotient fpone d))
(setq h (fptimes* h (fptimes* d c)))
(setq aa (fpquotient
(fpminus
(fptimes*
(fpplus a (intofp m))
(fptimes* (fpplus ab (intofp m)) z)))
(fptimes*
(fpplus a m2) (fpplus ap m2))))
(setq d (fpplus fpone (fptimes* aa d)))
(when (fplessp (fpabs d) beta-min) (setq d beta-min))
(setq c (fpplus fpone (fpquotient aa c)))
(when (fplessp (fpabs c) beta-min) (setq c beta-min))
(setq d (fpquotient fpone d))
(setq de (fptimes* d c))
(setq h (fptimes* h de))
(when (fplessp (fpabs (fpdifference de fpone)) beta-eps)
(return
(fpquotient
(fptimes*
h
(fptimes*
(cdr (exptbigfloat (bcons z) (bcons a)))
(cdr
(exptbigfloat (bcons (fpdifference fpone z)) (bcons b)))))
a)))))))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;