[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. d80256d66d585a69813f277b115d039cc620e62d
Raymond Toy
rtoy at common-lisp.net
Sat Mar 24 01:58:47 UTC 2012
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "OCT: A portable Lisp implementation for quad-double precision floats".
The branch, master has been updated
via d80256d66d585a69813f277b115d039cc620e62d (commit)
via b1ae6953934670284060ffc2810e424a5e0aac71 (commit)
via 82b854634c4531cab18173a265113c16e47fe1b3 (commit)
via 62f2b6b2074a593868d1b6fd0b42d7f85d3d5511 (commit)
via 19dd5231c42c6eec40a7876bd9ffb2661fb05012 (commit)
via 321ced44416bc5141b5306bc18f440eb179c4327 (commit)
via efbf11a6f5cd7e15f51cdb52fefd1584ab58c3d9 (commit)
via bc851c79e6fffe06c4383142bc2dfeb154fbde14 (commit)
via 6162b30bff1a78bcf7757476e8b05346af85f0e2 (commit)
via 7eb8fdd72b7227b93f95a0a90498b45a1caece0a (commit)
via 38b3e36708704f0de1bbfdc62bc6e52be08ab628 (commit)
via cc75066d522e0d1509f4a738908b913131c6deba (commit)
via f77890744f8e35f9c38ba19d435c6541f92d020b (commit)
via 06e250ea1734f4501d9c93289684c4bb5ae3c8f1 (commit)
via 3abcf3f7de7d89ba0fab0cb8f20931ebbfc6f844 (commit)
via dd77bf83d2f011da240e95f7eac29c97c9f0ef8d (commit)
from fe8cffb5ee8e7161addf7586ecaee00682c6bf1b (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit d80256d66d585a69813f277b115d039cc620e62d
Author: Raymond Toy <rtoy at google.com>
Date: Fri Mar 23 13:00:24 2012 -0700
Need to apply contagion for exp-integral-e.
diff --git a/qd-gamma.lisp b/qd-gamma.lisp
index 939b456..b233fab 100644
--- a/qd-gamma.lisp
+++ b/qd-gamma.lisp
@@ -528,27 +528,30 @@
;; for |arg(z)| < pi.
;;
;;
- (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))))
- ((< (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 (- 1 v) z)))
- (t
- ;; Use continued fraction for everything else.
- (cf-exp-integral-e v z))))
+ (incomplete-gamma-tail (- 1 v) z)))
+ (t
+ ;; Use continued fraction for everything else.
+ (cf-exp-integral-e v z)))))
;; Series for Fresnel S
;;
commit b1ae6953934670284060ffc2810e424a5e0aac71
Author: Raymond Toy <rtoy at google.com>
Date: Fri Mar 23 12:48:54 2012 -0700
Use exp-integral-e to evaluate incomplete-gamma-tail for real,
negative values of the parameter.
diff --git a/qd-gamma.lisp b/qd-gamma.lisp
index 6d80b71..939b456 100644
--- a/qd-gamma.lisp
+++ b/qd-gamma.lisp
@@ -375,9 +375,10 @@
(let* ((prec (float-contagion a z))
(a (apply-contagion a prec))
(z (apply-contagion z prec)))
- (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)))
;; For real values, we split the result to compute either the
commit 82b854634c4531cab18173a265113c16e47fe1b3
Author: Raymond Toy <rtoy at google.com>
Date: Fri Mar 23 10:38:31 2012 -0700
Document CHECK-ACCURACY
diff --git a/rt-tests.lisp b/rt-tests.lisp
index 8e116f8..d526af7 100644
--- a/rt-tests.lisp
+++ b/rt-tests.lisp
@@ -45,6 +45,12 @@
t
(- (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 (not (eq bits t))
commit 62f2b6b2074a593868d1b6fd0b42d7f85d3d5511
Author: Raymond Toy <rtoy at google.com>
Date: Fri Mar 23 10:21:40 2012 -0700
Bug fix for check-accuracy, and more tests for exp-integral-e.
* Don't let NaN's fool check-accuracy
* Add tests for exp-integral-e with v = 1.
diff --git a/rt-tests.lisp b/rt-tests.lisp
index 173c528..8e116f8 100644
--- a/rt-tests.lisp
+++ b/rt-tests.lisp
@@ -47,8 +47,10 @@
(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 est))
+ (not (float-nan-p bits))
+ (< bits limit))
(list bits limit est true)))))
(defvar *null* (make-broadcast-stream))
@@ -1477,4 +1479,19 @@
(e (exp-integral-e #q2 x))
(true #q0.326643862324553017730401565333637835828494690329010198058745549181386569998611289568))
(check-accuracy 208.4 e true))
- nil)
\ No newline at end of file
+ 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)
+
commit 19dd5231c42c6eec40a7876bd9ffb2661fb05012
Author: Raymond Toy <rtoy at google.com>
Date: Fri Mar 23 10:20:27 2012 -0700
Cleanups for psi and bug fix for exp-integral-e.
* Allow (exp-integral-e 1 z) to work.
* psi
* Handle psi(1) specially.
* Do a better job with cot(%pi*z) when z is an odd multiple of 1/2
where cot is 0.
* Fib bug in computing the number of terms when we try to float a
complex. Just float the realpart.
diff --git a/qd-gamma.lisp b/qd-gamma.lisp
index e353ff5..6d80b71 100644
--- a/qd-gamma.lisp
+++ b/qd-gamma.lisp
@@ -504,10 +504,11 @@
(- (psi v) (log z)))
(loop for k from 0
for term = 1 then (* term (/ -z k))
- for sum = (/ (- 1 v)) then (+ sum (let ((denom (- k n-1)))
- (if (zerop denom)
- 0
- (/ term denom))))
+ for sum = (if (= v 1) 0 (/ (- 1 v)))
+ then (+ sum (let ((denom (- k n-1)))
+ (if (zerop denom)
+ 0
+ (/ term denom))))
when (< (abs term) (* (abs sum) eps))
return sum)))
(loop for k from 0
@@ -798,18 +799,22 @@
;; So use reflection formula if Re(z) < 0. For z > 0, use the recurrence
;; 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 (- 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))
(x 0))
commit 321ced44416bc5141b5306bc18f440eb179c4327
Author: Raymond Toy <rtoy at google.com>
Date: Fri Mar 23 10:06:28 2012 -0700
Oops. In FLOAT, if it's already a float, don't change it if the
second arg is not given.
diff --git a/qd-methods.lisp b/qd-methods.lisp
index 10984a2..4fac522 100644
--- a/qd-methods.lisp
+++ b/qd-methods.lisp
@@ -289,7 +289,11 @@
(declaim (inline float))
(defun float (x &optional num-type)
- (qfloat x (or num-type 0.0)))
+ (if num-type
+ (qfloat x num-type)
+ (if (or (cl:floatp x) (typep x 'qd-real))
+ x
+ (qfloat x 0.0))))
(defmethod qrealpart ((x number))
(cl:realpart x))
commit efbf11a6f5cd7e15f51cdb52fefd1584ab58c3d9
Author: Raymond Toy <rtoy at google.com>
Date: Fri Mar 23 09:47:16 2012 -0700
Fix bug in FLOAT: second arg is optional! Add FLOAT-NAN-P method.
qd-methods.lisp:
* Second arg to {{{FLOAT}}} is optional.
* Add {{{FLOAT-NAN-P}}}.
qd-package.lisp:
* Need to shadow {{{EXT:FLOAT-NAN-P}}} on cmucl.
diff --git a/qd-methods.lisp b/qd-methods.lisp
index 0eef82d..10984a2 100644
--- a/qd-methods.lisp
+++ b/qd-methods.lisp
@@ -288,8 +288,8 @@
x)
(declaim (inline float))
-(defun float (x num-type)
- (qfloat x num-type))
+(defun float (x &optional num-type)
+ (qfloat x (or num-type 0.0)))
(defmethod qrealpart ((x number))
(cl:realpart x))
@@ -1111,6 +1111,13 @@ the same precision as the argument. The argument can be complex."))
(defmethod float-pi ((z qd-complex))
+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)
((function-name :accessor condition-function-name
diff --git a/qd-package.lisp b/qd-package.lisp
index 7db5cd3..7ca738b 100644
--- a/qd-package.lisp
+++ b/qd-package.lisp
@@ -211,6 +211,8 @@
#:rational
#:rationalize
)
+ #+cmu
+ (:shadow ext:float-nan-p)
;; Export types
(:export #:qd-real
#:qd-complex)
commit bc851c79e6fffe06c4383142bc2dfeb154fbde14
Merge: 6162b30 fe8cffb
Author: Raymond Toy <rtoy at google.com>
Date: Fri Mar 23 08:38:53 2012 -0700
Merge branch 'master' of git://common-lisp.net/projects/oct/oct
commit 6162b30bff1a78bcf7757476e8b05346af85f0e2
Merge: 7eb8fdd 4b332ed
Author: Raymond Toy <rtoy at google.com>
Date: Thu Mar 22 09:02:16 2012 -0700
Merge branch 'master' of git://common-lisp.net/projects/oct/oct
Conflicts:
qd-gamma.lisp
commit 7eb8fdd72b7227b93f95a0a90498b45a1caece0a
Author: Raymond Toy <rtoy at google.com>
Date: Wed Mar 21 16:20:04 2012 -0700
Allow exp-integral-e to work with negative integer orders. First
cut....
diff --git a/qd-gamma.lisp b/qd-gamma.lisp
index 9d87e1a..26a45a3 100644
--- a/qd-gamma.lisp
+++ b/qd-gamma.lisp
@@ -494,12 +494,28 @@
(let ((-z (- z))
(-v (- v))
(eps (epsilon z)))
- (loop for k from 0
- for term = 1 then (* term (/ -z k))
- for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
- when (< (abs term) (* (abs sum) eps))
- return (- (* (gamma (- 1 v)) (expt z (- v 1)))
- sum))))
+ (if (and (realp v)
+ (= v (ftruncate v)))
+ ;; v is an integer
+ (let ((n (truncate v)))
+ (- (* (/ (expt -z (- v 1))
+ (gamma v))
+ (- (psi v) (log z)))
+ (loop for k from 0 below n
+ for term = 1 then (* term (/ -z k))
+ for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
+ finally (return sum))
+ (loop for k from n
+ for term = 1 then (* term (/ -z k))
+ for sum = 0 then (+ sum (/ term (+ k 1 -v)))
+ when (< (abs term) (* (abs sum) eps))
+ return sum)))
+ (loop for k from 0
+ for term = 1 then (* term (/ -z k))
+ for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
+ when (< (abs term) (* (abs sum) eps))
+ return (- (* (gamma (- 1 v)) (expt z (- v 1)))
+ sum)))))
(defun exp-integral-e (v z)
"Exponential integral E:
commit 38b3e36708704f0de1bbfdc62bc6e52be08ab628
Author: Raymond Toy <rtoy at google.com>
Date: Wed Mar 21 15:29:07 2012 -0700
First cut at psi function.
diff --git a/qd-gamma.lisp b/qd-gamma.lisp
index 48f8345..9d87e1a 100644
--- a/qd-gamma.lisp
+++ b/qd-gamma.lisp
@@ -479,6 +479,16 @@
(* (- k)
(+ k v -1)))))))
+
+;; For v not an integer:
+;;
+;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf)
+;;
+;; For v an integer:
+;;
+;; E(v,z) = (-z)^(v-1)/(v-1)!*(psi(v)-log(z))
+;; - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf, k != n-1)
+;;
(defun s-exp-integral-e (v z)
;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf)
(let ((-z (- z))
@@ -679,3 +689,109 @@
(if (and (realp z) (plusp z))
(realpart (ci z))
(ci z))))
+
+(defconstant bern-values
+ (make-array 55
+ :initial-contents
+ '(1
+ -1/2
+ 1/6
+ 0
+ -1/30
+ 0
+ 1/42
+ 0
+ -1/30
+ 0
+ 5/66
+ 0
+ -691/2730
+ 0
+ 7/6
+ 0
+ -3617/510
+ 0
+ 43867/798
+ 0
+ -174611/330
+ 0
+ 854513/138
+ 0
+ -236364091/2730
+ 0
+ 8553103/6
+ 0
+ -23749461029/870
+ 0
+ 8615841276005/14322
+ 0
+ -7709321041217/510
+ 0
+ 2577687858367/6
+ 0
+ -26315271553053477373/1919190
+ 0
+ 2929993913841559/6
+ 0
+ -261082718496449122051/13530
+ 0
+ 1520097643918070802691/1806
+ 0
+ -27833269579301024235023/690
+ 0
+ 596451111593912163277961/282
+ 0
+ -5609403368997817686249127547/46410
+ 0
+ 495057205241079648212477525/66
+ 0
+ -801165718135489957347924991853/1590
+ 0
+ 29149963634884862421418123812691/798
+ )))
+
+(defun bern (k)
+ (aref bern-values k))
+
+(defun psi (z)
+ "Digamma function defined by
+
+ - %gamma + sum(1/k-1/(k+z-1), k, 1, inf)
+
+ where %gamma is Euler's constant"
+
+ ;; A&S 6.3.7: Reflection formula
+ ;;
+ ;; psi(1-z) = psi(z) + %pi*cot(%pi*z)
+ ;;
+ ;; A&S 6.3.6: Recurrence formula
+ ;;
+ ;; psi(n+z) = 1/(z+n-1)+1/(z+n-2)+...+1/(z+2)+1/(1+z)+psi(1+z)
+ ;;
+ ;; A&S 6.3.8: Asymptotic formula
+ ;;
+ ;; psi(z) ~ log(z) - sum(bern(2*n)/(2*n*z^(2*n)), n, 1, inf)
+ ;;
+ ;; So use reflection formula if Re(z) < 0. For z > 0, use the recurrence
+ ;; formula to increase the argument and then apply the asymptotic formula.
+
+ (cond ((minusp (realpart z))
+ (- (psi (- 1 z))
+ (* +pi+ (/ (tan (* +pi+ z))))))
+ (t
+ (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon z) 10)))))))
+ (m 0)
+ (y (expt (+ z k) 2))
+ (x 0))
+ (loop for i from 1 upto (floor k 2) do
+ (progn
+ (incf m (+ (/ (+ z i i -1))
+ (/ (+ z i i -2))))
+ (setf x (/ (+ x (/ (bern (+ k 2 (* -2 i)))
+ (- k i i -2)))
+ y))))
+ (- (log (+ z k))
+ (/ (* 2 (+ z k)))
+ x
+ m)))))
+
commit cc75066d522e0d1509f4a738908b913131c6deba
Merge: f778907 405df61
Author: Raymond Toy <rtoy at google.com>
Date: Wed Mar 21 09:37:15 2012 -0700
Merge branch 'master' of git://common-lisp.net/projects/oct/oct
commit f77890744f8e35f9c38ba19d435c6541f92d020b
Merge: 06e250e c388f81
Author: Raymond Toy <rtoy at google.com>
Date: Tue Dec 6 08:50:45 2011 -0800
Merge branch 'master' of git://common-lisp.net/projects/oct/oct
commit 06e250ea1734f4501d9c93289684c4bb5ae3c8f1
Merge: 3abcf3f 26ef83a
Author: Raymond Toy <rtoy at google.com>
Date: Mon Dec 5 14:07:49 2011 -0800
Merge branch 'master' of git://common-lisp.net/projects/oct/oct
commit 3abcf3f7de7d89ba0fab0cb8f20931ebbfc6f844
Merge: dd77bf8 06f8a53
Author: Raymond Toy <rtoy at google.com>
Date: Mon Dec 5 13:44:51 2011 -0800
Merge branch 'master' of git://common-lisp.net/projects/oct/oct
commit dd77bf83d2f011da240e95f7eac29c97c9f0ef8d
Author: Raymond Toy <rtoy at google.com>
Date: Mon Dec 5 09:43:15 2011 -0800
Update asdf version (and make it compatible with latest asdf).
diff --git a/oct.asd b/oct.asd
index a6111ca..d6e7090 100644
--- a/oct.asd
+++ b/oct.asd
@@ -36,7 +36,7 @@
:author "Raymond Toy"
:maintainer "See <http://www.common-lisp.net/project/oct>"
:licence "MIT"
- :version "2011-02-09" ; Just use the date
+ :version "2011.12.05" ; Just use the date
:components
((:file "qd-package")
(:file "qd-rep" :depends-on ("qd-package"))
@@ -74,7 +74,7 @@
(asdf:defsystem oct-tests
:depends-on (oct)
- :version "2011-02-09" ; Just use the date
+ :version "2011.12.05" ; Just use the date
:in-order-to ((compile-op (load-op :rt))
(test-op (load-op :rt :oct)))
:components
-----------------------------------------------------------------------
Summary of changes:
qd-gamma.lisp | 79 ++++++++++++++++++++++++++++++------------------------
qd-methods.lisp | 15 +++++++++-
qd-package.lisp | 2 +
rt-tests.lisp | 29 ++++++++++++++++++--
4 files changed, 85 insertions(+), 40 deletions(-)
hooks/post-receive
--
OCT: A portable Lisp implementation for quad-double precision floats
More information about the oct-cvs
mailing list