Changeset 6cfb0a


Ignore:
Timestamp:
04/09/12 15:37:51 (3 years ago)
Author:
Raymond Toy <toy.raymond@…>
Branches:
master
Children:
cb1a5d
Parents:
8ec0d2 (diff), dd5c2a (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

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

Files:
1 added
5 edited

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  
    11331133  +pi+)
    11341134
     1135(defmethod float-nan-p ((x cl:float))
     1136  ;; CMUCL has ext:float-nan-p.  Should we use that instead?
     1137  (not (= x x)))
     1138
     1139(defmethod float-nan-p ((x qd-real))
     1140  (float-nan-p (qd-parts (qd-value x))))
     1141
    11351142
    11361143
  • 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
     
    15001508      (check-accuracy 208.4 e true))
    15011509  nil)
     1510
     1511(rt:deftest expintegral-e.6d
     1512    (let* ((x .5d0)
     1513           (e (exp-integral-e 1d0 x))
     1514           (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883))
     1515      (check-accuracy 53.9 e true))
     1516  nil)
     1517
     1518(rt:deftest expintegral-e.6q
     1519    (let* ((x #q.5)
     1520           (e (exp-integral-e #q1 x))
     1521           (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883))
     1522      (check-accuracy 219.1 e true))
     1523  nil)
     1524
Note: See TracChangeset for help on using the changeset viewer.