Changeset 6cfb0a
 Timestamp:
 04/09/12 15:37:51 (3 years ago)
 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.  Files:

 1 added
 5 edited
Legend:
 Unmodified
 Added
 Removed

oct.asd
r405df6 rdbc1e3 68 68 (:file "qdgamma" 69 69 :dependson ("qdmethods" "qdreader")) 70 )) 70 (:file "qdbessel" 71 :dependson ("qdmethods")))) 71 72 72 73 (defmethod perform ((op testop) (c (eql (asdf:findsystem :oct)))) 
qdgamma.lisp
r6cd962 r6cfb0a 380 380 integrate(t^(a1)*exp(t), t, z, inf)" 381 381 (withfloatingpointcontagion (a z) 382 (if (zerop a) 383 ;; incomplete_gamma_tail(0, z) = exp_integral_e(1,z) 384 (expintegrale 1 z) 382 (if (and (realp a) (<= a 0)) 383 ;; incomplete_gamma_tail(v, z) = z^v*exp_integral_e(1a,z) 384 (* (expt z a) 385 (expintegrale ( 1 a) z)) 385 386 (if (and (zerop (imagpart a)) 386 387 (zerop (imagpart z))) … … 532 533 ;; 533 534 ;; 534 (cond ((and (realp v) (minusp v)) 535 ;; E(v, z) = z^(v1)*incomplete_gamma_tail(v+1,z) 536 (let ((v ( v))) 535 (let* ((prec (floatcontagion v z)) 536 (v (applycontagion v prec)) 537 (z (applycontagion z prec))) 538 (cond ((and (realp v) (minusp v)) 539 ;; E(v, z) = z^(v1)*incomplete_gamma_tail(v+1,z) 540 (let ((v ( v))) 541 (* (expt z ( v 1)) 542 (incompletegammatail (+ v 1) z)))) 543 ((< (abs z) 1) 544 ;; Use series for small z 545 (sexpintegrale 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 incompletegammatail 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^(v1)*incomplete_gamma_tail(1v,z) 537 554 (* (expt z ( v 1)) 538 555 (incompletegammatail (+ v 1) z)))) … … 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 (floatpi z))) 801 821 (flet ((cotpi (z) 802 ;; cot(%pi*z), car 803 (handlercase 804 (/ (tan (* p z))) 805 (divisionbyzero () 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 (cotpi 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)) 
qdmethods.lisp
r8ec0d2 r6cfb0a 1133 1133 +pi+) 1134 1134 1135 (defmethod floatnanp ((x cl:float)) 1136 ;; CMUCL has ext:floatnanp. Should we use that instead? 1137 (not (= x x))) 1138 1139 (defmethod floatnanp ((x qdreal)) 1140 (floatnanp (qdparts (qdvalue x)))) 1141 1135 1142 1136 1143 
qdpackage.lisp
r8ec0d2 r6cfb0a 213 213 #:rationalize 214 214 ) 215 #+cmu 216 (:shadow ext:floatnanp) 215 217 ;; Export types 216 218 (:export #:qdreal 
rttests.lisp
r8ec0d2 r6cfb0a 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 checkaccuracy (limit est true) 49 55 (let ((bits (bitaccuracy est true))) 50 (if (numberp bits) 51 (if (< bits limit) 56 (if (not (eq bits t)) 57 (if (and (not (floatnanp (realpart est))) 58 (not (floatnanp bits)) 59 (< bits limit)) 52 60 (list bits limit est true))))) 53 61 … … 1500 1508 (checkaccuracy 208.4 e true)) 1501 1509 nil) 1510 1511 (rt:deftest expintegrale.6d 1512 (let* ((x .5d0) 1513 (e (expintegrale 1d0 x)) 1514 (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883)) 1515 (checkaccuracy 53.9 e true)) 1516 nil) 1517 1518 (rt:deftest expintegrale.6q 1519 (let* ((x #q.5) 1520 (e (expintegrale #q1 x)) 1521 (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883)) 1522 (checkaccuracy 219.1 e true)) 1523 nil) 1524
Note: See TracChangeset
for help on using the changeset viewer.