Changeset 6cfb0ac4b6bcc1a25bc119e87fd2b57bfa1f4355

Timestamp:

04/09/12 08:37:51 (14 months ago)

Author:

Raymond Toy <toy.raymond@…>

Children:

cb1a5d41baf9d4db12a2563a230d0a1e55c8adea

Parents:

8ec0d2004ed4063fd568250b00d940b208573a92, dd5c2a4a6628648c5fe9a7eec98c82a9d284daa3

git-committer:

Raymond Toy <toy.raymond@…> (04/09/12 08:37:51)

Message:

Merge branch 'master' of  ssh://common-lisp.net/var/git/projects/oct/oct

Files:

• oct.asd (modified) (1 diff)
• qd-bessel.lisp (added)
• qd-gamma.lisp (modified) (3 diffs)
• qd-methods.lisp (modified) (1 diff)
• qd-package.lisp (modified) (1 diff)
• rt-tests.lisp (modified) (2 diffs)

oct.asd

@@ -68 +68 @@
    (:file "qd-gamma"
           :depends-on ("qd-methods" "qd-reader"))
-   ))
+   (:file "qd-bessel"
+          :depends-on ("qd-methods"))))
 
 (defmethod perform ((op test-op) (c (eql (asdf:find-system :oct))))

qd-gamma.lisp

@@ -380 +380 @@
   integrate(t^(a-1)*exp(-t), t, z, inf)"
   (with-floating-point-contagion (a z)
-    (if (zerop a)
-        ;; incomplete_gamma_tail(0, z) = exp_integral_e(1,z)
-        (exp-integral-e 1 z)
+    (if (and (realp a) (<= a 0))
+        ;; incomplete_gamma_tail(v, z) = z^v*exp_integral_e(1-a,z)
+        (* (expt z a)
+           (exp-integral-e (- 1 a) z))
         (if (and (zerop (imagpart a))
                  (zerop (imagpart z)))
@@ -532 +533 @@
   ;;
   ;;
-  (cond ((and (realp v) (minusp v))
-         ;; E(-v, z) = z^(-v-1)*incomplete_gamma_tail(v+1,z)
-         (let ((-v (- v)))
+  (let* ((prec (float-contagion v z))
+         (v (apply-contagion v prec))
+         (z (apply-contagion z prec)))
+    (cond ((and (realp v) (minusp v))
+           ;; E(-v, z) = z^(-v-1)*incomplete_gamma_tail(v+1,z)
+           (let ((-v (- v)))
+             (* (expt z (- v 1))
+                (incomplete-gamma-tail (+ -v 1) z))))
+          ((< (abs z) 1)
+           ;; Use series for small z
+           (s-exp-integral-e v z))
+          ((>= (abs (phase z)) 3.1)
+           ;; The continued fraction doesn't converge on the negative
+           ;; real axis, and converges very slowly near the negative
+           ;; real axis, so use the incomplete-gamma-tail function in
+           ;; this region.  "Closeness" to the negative real axis is
+           ;; teken to mean that z is in a sector near the axis.
+           ;;
+           ;; E(v,z) = z^(v-1)*incomplete_gamma_tail(1-v,z)
            (* (expt z (- v 1))
               (incomplete-gamma-tail (+ -v 1) z))))
@@ -797 +814 @@
   ;; formula to increase the argument and then apply the asymptotic formula.
 
-  (cond ((minusp (realpart z))
-         (let ((p (float +pi+ (realpart z))))
+  (cond ((= z 1)
+         ;; psi(1) = -%gamma
+         (- (float +%gamma+ (if (integerp z) 0.0 z))))
+        ((minusp (realpart z))
+         (let ((p (float-pi z)))
            (flet ((cot-pi (z)
-                    ;; cot(%pi*z), car
-                    (handler-case
-                        (/ (tan (* p z)))
-                      (division-by-zero ()
-                        (* 0 z)))))
+                    ;; cot(%pi*z), carefully.  If z is an odd multiple
+                    ;; of 1/2, cot is 0.
+                    (if (and (realp z)
+                             (= 1/2 (abs (- z (ftruncate z)))))
+                        (float 0 z)
+                        (/ (tan (* p z))))))
              (- (psi (- 1 z))
                 (* p (cot-pi z))))))
         (t
-         (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon z) 10)))))))
+         (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon (float (realpart z))) 10)))))))
                 (m 0)
                 (y (expt (+ z k) 2))

qd-methods.lisp

@@ -1130 +1130 @@
   +pi+)
 
+(defmethod float-nan-p ((x cl:float))
+  ;; CMUCL has ext:float-nan-p.  Should we use that instead?
+  (not (= x x)))
+
+(defmethod float-nan-p ((x qd-real))
+  (float-nan-p (qd-parts (qd-value x))))
+
 
 (define-condition domain-error (simple-error)

qd-package.lisp

@@ -213 +213 @@
            #:rationalize
            )
+  #+cmu
+  (:shadow ext:float-nan-p)
   ;; Export types
   (:export #:qd-real

rt-tests.lisp

@@ -46 +46 @@
         (- (log err 2)))))
 
+;; Check actual value EST is with LIMIT bits of the true value TRUE.
+;; If so, return NIL.  Otherwise, return a list of the actual bits of
+;; accuracy, the desired accuracy, and the values.  This is mostly to
+;; make it easy to see what the actual accuracy was and the arguments
+;; for the test, which is important for the tests that use random
+;; values.
 (defun check-accuracy (limit est true)
   (let ((bits (bit-accuracy est true)))
-    (if (numberp bits)
-        (if (< bits limit)
+    (if (not (eq bits t))
+        (if (and (not (float-nan-p (realpart est)))
+                 (not (float-nan-p bits))
+                 (< bits limit))
             (list bits limit est true)))))
 
@@ -1495 +1503 @@
       (check-accuracy 208.4 e true))
   nil)
+
+(rt:deftest expintegral-e.6d
+    (let* ((x .5d0)
+           (e (exp-integral-e 1d0 x))
+           (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883))
+      (check-accuracy 53.9 e true))
+  nil)
+
+(rt:deftest expintegral-e.6q
+    (let* ((x #q.5)
+           (e (exp-integral-e #q1 x))
+           (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883))
+      (check-accuracy 219.1 e true))
+  nil)
+