Changeset 6cd962


Ignore:
Timestamp:
04/08/12 16:35:35 (3 years ago)
Author:
Raymond Toy <toy.raymond@…>
Branches:
master
Children:
dbc1e3
Parents:
c86217
Message:

Define macro WITH-FLOATING-POINT-CONTAGION.

qd-methods.lisp:
o Define the macro

qd-gamma.lisp:
o Use it.

Files:
2 edited

Legend:

Unmodified
Added
Removed
  • qd-gamma.lisp

    rc86217 r6cd962  
    109109           ;; Or
    110110           ;;   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)))))
    115116          (t
    116117           (let ((absz (abs z)))
     
    378379
    379380  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)
    383382    (if (zerop a)
    384383        ;; incomplete_gamma_tail(0, z) = exp_integral_e(1,z)
     
    410409
    411410  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)
    415412    (if (and (< (abs a) 1) (< (abs z) 1))
    416413        (s-incomplete-gamma a z)
     
    540537           (* (expt z (- v 1))
    541538              (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.
    544544         (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)))
    555545        (t
    556546         ;; Use continued fraction for everything else.
  • qd-methods.lisp

    r06f8a5 r6cd962  
    8080              (coerce (imagpart number) precision)))))
    8181
     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))))
    8293
    8394(defmethod add1 ((a number))
Note: See TracChangeset for help on using the changeset viewer.