Changeset 6cfb0ac4b6bcc1a25bc119e87fd2b57bfa1f4355

Show
Ignore:
Timestamp:
04/09/12 08:37:51 (2 years 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:
1 added
5 modified

Legend:

Unmodified
Added
Removed
  • oct.asd

    r405df6 rdbc1e3  
    6868   (:file "qd-gamma" 
    6969          :depends-on ("qd-methods" "qd-reader")) 
    70    )) 
     70   (:file "qd-bessel" 
     71          :depends-on ("qd-methods")))) 
    7172 
    7273(defmethod perform ((op test-op) (c (eql (asdf:find-system :oct)))) 
  • qd-gamma.lisp

    r6cd962 r6cfb0a  
    380380  integrate(t^(a-1)*exp(-t), t, z, inf)" 
    381381  (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)) 
    385386        (if (and (zerop (imagpart a)) 
    386387                 (zerop (imagpart z))) 
     
    532533  ;; 
    533534  ;; 
    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) 
    537554           (* (expt z (- v 1)) 
    538555              (incomplete-gamma-tail (+ -v 1) z)))) 
     
    797814  ;; formula to increase the argument and then apply the asymptotic formula. 
    798815 
    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))) 
    801821           (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)))))) 
    807828             (- (psi (- 1 z)) 
    808829                (* p (cot-pi z)))))) 
    809830        (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))))))) 
    811832                (m 0) 
    812833                (y (expt (+ z k) 2)) 
  • qd-methods.lisp

    r8ec0d2 r6cfb0a  
    11301130  +pi+) 
    11311131 
     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 
    11321139 
    11331140(define-condition domain-error (simple-error) 
  • qd-package.lisp

    r8ec0d2 r6cfb0a  
    213213           #:rationalize 
    214214           ) 
     215  #+cmu 
     216  (:shadow ext:float-nan-p) 
    215217  ;; Export types 
    216218  (:export #:qd-real 
  • rt-tests.lisp

    r8ec0d2 r6cfb0a  
    4646        (- (log err 2))))) 
    4747 
     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. 
    4854(defun check-accuracy (limit est true) 
    4955  (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)) 
    5260            (list bits limit est true))))) 
    5361 
     
    14951503      (check-accuracy 208.4 e true)) 
    14961504  nil) 
     1505 
     1506(rt:deftest expintegral-e.6d 
     1507    (let* ((x .5d0) 
     1508           (e (exp-integral-e 1d0 x)) 
     1509           (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883)) 
     1510      (check-accuracy 53.9 e true)) 
     1511  nil) 
     1512 
     1513(rt:deftest expintegral-e.6q 
     1514    (let* ((x #q.5) 
     1515           (e (exp-integral-e #q1 x)) 
     1516           (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883)) 
     1517      (check-accuracy 219.1 e true)) 
     1518  nil) 
     1519