Commit 6cd96249b94dc29962cecda996a5799d8ac27271

Authored by Raymond Toy
1 parent c86217a5

Define macro WITH-FLOATING-POINT-CONTAGION.

qd-methods.lisp:
o Define the macro

qd-gamma.lisp:
o Use it.
Showing 2 changed files with 23 additions and 22 deletions
qd-gamma.lisp
... ... @@ -108,10 +108,11 @@
108 108 ;; log(gamma(-z)) = log(pi)-log(-z)-log(sin(pi*z))-log(gamma(z))
109 109 ;; Or
110 110 ;; log(gamma(z)) = log(pi)-log(-z)-log(sin(pi*z))-log(gamma(-z))
111   - (- (apply-contagion (log pi) precision)
112   - (log (- z))
113   - (apply-contagion (log (sin (* pi z))) precision)
114   - (log-gamma (- z))))
  111 + (let ((p (float-pi z)))
  112 + (- (log p)
  113 + (log (- z))
  114 + (log (sin (* p z)))
  115 + (log-gamma (- z)))))
115 116 (t
116 117 (let ((absz (abs z)))
117 118 (cond ((>= absz limit)
... ... @@ -377,9 +378,7 @@
377 378 "Tail of the incomplete gamma function defined by:
378 379  
379 380 integrate(t^(a-1)*exp(-t), t, z, inf)"
380   - (let* ((prec (float-contagion a z))
381   - (a (apply-contagion a prec))
382   - (z (apply-contagion z prec)))
  381 + (with-floating-point-contagion (a z)
383 382 (if (zerop a)
384 383 ;; incomplete_gamma_tail(0, z) = exp_integral_e(1,z)
385 384 (exp-integral-e 1 z)
... ... @@ -409,9 +408,7 @@
409 408 "Incomplete gamma function defined by:
410 409  
411 410 integrate(t^(a-1)*exp(-t), t, 0, z)"
412   - (let* ((prec (float-contagion a z))
413   - (a (apply-contagion a prec))
414   - (z (apply-contagion z prec)))
  411 + (with-floating-point-contagion (a z)
415 412 (if (and (< (abs a) 1) (< (abs z) 1))
416 413 (s-incomplete-gamma a z)
417 414 (if (and (realp a) (realp z))
... ... @@ -539,19 +536,12 @@
539 536 (let ((-v (- v)))
540 537 (* (expt z (- v 1))
541 538 (incomplete-gamma-tail (+ -v 1) z))))
542   - ((< (abs z) 1)
543   - ;; Use series for small z
  539 + ((or (< (abs z) 1) (>= (abs (phase z)) 3.1))
  540 + ;; Use series for small z or if z is near the negative real
  541 + ;; axis because the continued fraction does not converge on
  542 + ;; the negative axis and converges slowly near the negative
  543 + ;; axis.
544 544 (s-exp-integral-e v z))
545   - ((>= (abs (phase z)) 3.1)
546   - ;; The continued fraction doesn't converge on the negative
547   - ;; real axis, and converges very slowly near the negative
548   - ;; real axis, so use the incomplete-gamma-tail function in
549   - ;; this region. "Closeness" to the negative real axis is
550   - ;; teken to mean that z is in a sector near the axis.
551   - ;;
552   - ;; E(v,z) = z^(v-1)*incomplete_gamma_tail(1-v,z)
553   - (* (expt z (- v 1))
554   - (incomplete-gamma-tail (- 1 v) z)))
555 545 (t
556 546 ;; Use continued fraction for everything else.
557 547 (cf-exp-integral-e v z))))
... ...
qd-methods.lisp
... ... @@ -79,6 +79,17 @@
79 79 (complex (coerce (realpart number) precision)
80 80 (coerce (imagpart number) precision)))))
81 81  
  82 +;; WITH-FLOATING-POINT-CONTAGION - macro
  83 +;;
  84 +;; Determines the highest precision of the variables in VARLIST and
  85 +;; converts each of the values to that precision.
  86 +(defmacro with-floating-point-contagion (varlist &body body)
  87 + (let ((precision (gensym "PRECISION-")))
  88 + `(let ((,precision (float-contagion ,@varlist)))
  89 + (let (,@(mapcar #'(lambda (v)
  90 + `(,v (apply-contagion ,v ,precision)))
  91 + varlist))
  92 + ,@body))))
82 93  
83 94 (defmethod add1 ((a number))
84 95 (cl::1+ a))
... ...