Commit 6cfb0ac4b6bcc1a25bc119e87fd2b57bfa1f4355

Authored by Raymond Toy
2 parents 8ec0d200 dd5c2a4a
qd-gamma.lisp
... ... @@ -379,9 +379,10 @@
379 379  
380 380 integrate(t^(a-1)*exp(-t), t, z, inf)"
381 381 (with-floating-point-contagion (a z)
382   - (if (zerop a)
383   - ;; incomplete_gamma_tail(0, z) = exp_integral_e(1,z)
384   - (exp-integral-e 1 z)
  382 + (if (and (realp a) (<= a 0))
  383 + ;; incomplete_gamma_tail(v, z) = z^v*exp_integral_e(1-a,z)
  384 + (* (expt z a)
  385 + (exp-integral-e (- 1 a) z))
385 386 (if (and (zerop (imagpart a))
386 387 (zerop (imagpart z)))
387 388 ;; For real values, we split the result to compute either the
... ... @@ -531,9 +532,25 @@
531 532 ;; for |arg(z)| < pi.
532 533 ;;
533 534 ;;
534   - (cond ((and (realp v) (minusp v))
535   - ;; E(-v, z) = z^(-v-1)*incomplete_gamma_tail(v+1,z)
536   - (let ((-v (- v)))
  535 + (let* ((prec (float-contagion v z))
  536 + (v (apply-contagion v prec))
  537 + (z (apply-contagion z prec)))
  538 + (cond ((and (realp v) (minusp v))
  539 + ;; E(-v, z) = z^(-v-1)*incomplete_gamma_tail(v+1,z)
  540 + (let ((-v (- v)))
  541 + (* (expt z (- v 1))
  542 + (incomplete-gamma-tail (+ -v 1) z))))
  543 + ((< (abs z) 1)
  544 + ;; Use series for small z
  545 + (s-exp-integral-e v z))
  546 + ((>= (abs (phase z)) 3.1)
  547 + ;; The continued fraction doesn't converge on the negative
  548 + ;; real axis, and converges very slowly near the negative
  549 + ;; real axis, so use the incomplete-gamma-tail function in
  550 + ;; this region. "Closeness" to the negative real axis is
  551 + ;; teken to mean that z is in a sector near the axis.
  552 + ;;
  553 + ;; E(v,z) = z^(v-1)*incomplete_gamma_tail(1-v,z)
537 554 (* (expt z (- v 1))
538 555 (incomplete-gamma-tail (+ -v 1) z))))
539 556 ((or (< (abs z) 1) (>= (abs (phase z)) 3.1))
... ... @@ -796,18 +813,22 @@
796 813 ;; So use reflection formula if Re(z) < 0. For z > 0, use the recurrence
797 814 ;; formula to increase the argument and then apply the asymptotic formula.
798 815  
799   - (cond ((minusp (realpart z))
800   - (let ((p (float +pi+ (realpart z))))
  816 + (cond ((= z 1)
  817 + ;; psi(1) = -%gamma
  818 + (- (float +%gamma+ (if (integerp z) 0.0 z))))
  819 + ((minusp (realpart z))
  820 + (let ((p (float-pi z)))
801 821 (flet ((cot-pi (z)
802   - ;; cot(%pi*z), car
803   - (handler-case
804   - (/ (tan (* p z)))
805   - (division-by-zero ()
806   - (* 0 z)))))
  822 + ;; cot(%pi*z), carefully. If z is an odd multiple
  823 + ;; of 1/2, cot is 0.
  824 + (if (and (realp z)
  825 + (= 1/2 (abs (- z (ftruncate z)))))
  826 + (float 0 z)
  827 + (/ (tan (* p z))))))
807 828 (- (psi (- 1 z))
808 829 (* p (cot-pi z))))))
809 830 (t
810   - (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon z) 10)))))))
  831 + (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon (float (realpart z))) 10)))))))
811 832 (m 0)
812 833 (y (expt (+ z k) 2))
813 834 (x 0))
... ...
qd-methods.lisp
... ... @@ -1129,6 +1129,13 @@ the same precision as the argument. The argument can be complex.&quot;))
1129 1129 (defmethod float-pi ((z qd-complex))
1130 1130 +pi+)
1131 1131  
  1132 +(defmethod float-nan-p ((x cl:float))
  1133 + ;; CMUCL has ext:float-nan-p. Should we use that instead?
  1134 + (not (= x x)))
  1135 +
  1136 +(defmethod float-nan-p ((x qd-real))
  1137 + (float-nan-p (qd-parts (qd-value x))))
  1138 +
1132 1139  
1133 1140 (define-condition domain-error (simple-error)
1134 1141 ((function-name :accessor condition-function-name
... ... @@ -1139,4 +1146,4 @@ the same precision as the argument. The argument can be complex.&quot;))
1139 1146 (pprint-logical-block (stream nil :per-line-prefix " ")
1140 1147 (apply #'format stream
1141 1148 (simple-condition-format-control condition)
1142   - (simple-condition-format-arguments condition))))))
1143 1149 \ No newline at end of file
  1150 + (simple-condition-format-arguments condition))))))
... ...
qd-package.lisp
... ... @@ -212,6 +212,8 @@
212 212 #:rational
213 213 #:rationalize
214 214 )
  215 + #+cmu
  216 + (:shadow ext:float-nan-p)
215 217 ;; Export types
216 218 (:export #:qd-real
217 219 #:qd-complex)
... ...
rt-tests.lisp
... ... @@ -45,10 +45,18 @@
45 45 t
46 46 (- (log err 2)))))
47 47  
  48 +;; Check actual value EST is with LIMIT bits of the true value TRUE.
  49 +;; If so, return NIL. Otherwise, return a list of the actual bits of
  50 +;; accuracy, the desired accuracy, and the values. This is mostly to
  51 +;; make it easy to see what the actual accuracy was and the arguments
  52 +;; for the test, which is important for the tests that use random
  53 +;; values.
48 54 (defun check-accuracy (limit est true)
49 55 (let ((bits (bit-accuracy est true)))
50   - (if (numberp bits)
51   - (if (< bits limit)
  56 + (if (not (eq bits t))
  57 + (if (and (not (float-nan-p (realpart est)))
  58 + (not (float-nan-p bits))
  59 + (< bits limit))
52 60 (list bits limit est true)))))
53 61  
54 62 (defvar *null* (make-broadcast-stream))
... ... @@ -1493,4 +1501,19 @@
1493 1501 (e (exp-integral-e #q2 x))
1494 1502 (true #q0.326643862324553017730401565333637835828494690329010198058745549181386569998611289568))
1495 1503 (check-accuracy 208.4 e true))
1496   - nil)
1497 1504 \ No newline at end of file
  1505 + nil)
  1506 +
  1507 +(rt:deftest expintegral-e.6d
  1508 + (let* ((x .5d0)
  1509 + (e (exp-integral-e 1d0 x))
  1510 + (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883))
  1511 + (check-accuracy 53.9 e true))
  1512 + nil)
  1513 +
  1514 +(rt:deftest expintegral-e.6q
  1515 + (let* ((x #q.5)
  1516 + (e (exp-integral-e #q1 x))
  1517 + (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883))
  1518 + (check-accuracy 219.1 e true))
  1519 + nil)
  1520 +
... ...