Changeset 6cd96249b94dc29962cecda996a5799d8ac27271

Show
Ignore:
Timestamp:
04/08/12 09:35:35 (2 years ago)
Author:
Raymond Toy <toy.raymond@…>
Children:
dbc1e376add1d492dc35c37b3098c2ffb796d6f1
Parents:
c86217a5f7191f1a29566d6ee24cb57344dd546d
git-committer:
Raymond Toy <toy.raymond@…> (04/08/12 09:35:35)
Message:

Define macro WITH-FLOATING-POINT-CONTAGION.

qd-methods.lisp:
o Define the macro

qd-gamma.lisp:
o Use it.

Files:
2 modified

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))