March 2012 - oct-cvs - mailman3.common-lisp.net

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. dd5c2a4a6628648c5fe9a7eec98c82a9d284daa3
by Raymond Toy 24 Mar '12

24 Mar '12

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 dd5c2a4a6628648c5fe9a7eec98c82a9d284daa3 (commit) via 41f1ca3edadee2d6e00ba7dc287ea66d92897ca6 (commit) from d80256d66d585a69813f277b115d039cc620e62d (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 dd5c2a4a6628648c5fe9a7eec98c82a9d284daa3 Author: Raymond Toy <rtoy(a)google.com> Date: Fri Mar 23 19:10:16 2012 -0700 Need to take ABS before testing for 1/2. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index b233fab..f0e2418 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -812,7 +812,7 @@ ;; 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)))) + (= 1/2 (abs (- z (ftruncate z))))) (float 0 z) (/ (tan (* p z)))))) (- (psi (- 1 z)) commit 41f1ca3edadee2d6e00ba7dc287ea66d92897ca6 Author: Raymond Toy <rtoy(a)google.com> Date: Fri Mar 23 19:05:24 2012 -0700 Oops. Need to take the realpart before call FLOAT-NAN-P. diff --git a/rt-tests.lisp b/rt-tests.lisp index d526af7..8732e28 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -54,7 +54,7 @@ (defun check-accuracy (limit est true) (let ((bits (bit-accuracy est true))) (if (not (eq bits t)) - (if (and (not (float-nan-p est)) + (if (and (not (float-nan-p (realpart est))) (not (float-nan-p bits)) (< bits limit)) (list bits limit est true))))) ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 2 +- rt-tests.lisp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. d80256d66d585a69813f277b115d039cc620e62d
by Raymond Toy 24 Mar '12

24 Mar '12

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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. fe8cffb5ee8e7161addf7586ecaee00682c6bf1b
by Raymond Toy 23 Mar '12

23 Mar '12

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 fe8cffb5ee8e7161addf7586ecaee00682c6bf1b (commit) from e10e402d446339130042607357b0129b35a1faa7 (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 fe8cffb5ee8e7161addf7586ecaee00682c6bf1b Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Thu Mar 22 22:58:43 2012 -0700 Clean up implemenation of s-exp-integral-e a bit. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index 13f0bd9..e353ff5 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -497,18 +497,19 @@ (if (and (realp v) (= v (ftruncate v))) ;; v is an integer - (let ((n (truncate v))) - (- (* (/ (expt -z (- v 1)) + (let* ((n (truncate v)) + (n-1 (1- n))) + (- (* (/ (expt -z n-1) (gamma v)) (- (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 1 (- n)))) + for term = 1 then (* term (/ -z k)) + for sum = (/ (- 1 v)) then (+ sum (let ((denom (- k n-1))) (if (zerop denom) 0 - (/ term (+ k 1 -v))))) - when (< (abs term) (* (abs sum) eps)) - return sum))) + (/ term denom)))) + 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))) ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 15 ++++++++------- 1 files changed, 8 insertions(+), 7 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. e10e402d446339130042607357b0129b35a1faa7
by Raymond Toy 23 Mar '12

23 Mar '12

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 e10e402d446339130042607357b0129b35a1faa7 (commit) via 1d9ec007bb4172eaccb8ca1db543c1218192cdb5 (commit) from 4b332ed2140e56c6fcaa689e6b8a48be36c988f6 (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 e10e402d446339130042607357b0129b35a1faa7 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Thu Mar 22 21:13:07 2012 -0700 Correct the computation of s-exp-integral-e when v is an integer. Add tests. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index e5fc5f0..13f0bd9 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -501,16 +501,14 @@ (- (* (/ (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))) - when (< (abs term) (* (abs sum) eps)) - 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 (let ((denom (+ k 1 (- n)))) + (if (zerop denom) + 0 + (/ 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))) diff --git a/rt-tests.lisp b/rt-tests.lisp index fa6da81..173c528 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1464,3 +1464,17 @@ (expt x 4)))) (check-accuracy 210.8 e true)) nil) + +(rt:deftest expintegral-e.5d + (let* ((x .5d0) + (e (exp-integral-e 2d0 x)) + (true #q0.3266438623245530177304015653336378358284946903290101)) + (check-accuracy 51.2 e true)) + nil) + +(rt:deftest expintegral-e.5q + (let* ((x #q.5) + (e (exp-integral-e #q2 x)) + (true #q0.326643862324553017730401565333637835828494690329010198058745549181386569998611289568)) + (check-accuracy 208.4 e true)) + nil) \ No newline at end of file commit 1d9ec007bb4172eaccb8ca1db543c1218192cdb5 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Wed Mar 21 22:24:48 2012 -0700 Handle the case of negative v in exp-integral-e(v, z). Add tests. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index d49acf8..e5fc5f0 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -527,7 +527,12 @@ ;; for |arg(z)| < pi. ;; ;; - (cond ((< (abs z) 1) + (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) diff --git a/rt-tests.lisp b/rt-tests.lisp index c5b6743..fa6da81 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1447,4 +1447,20 @@ (let* ((e (exp-integral-e 2 1d0)) (true 0.14849550677592204791835999d0)) (check-accuracy 47.5 e true)) - nil) \ No newline at end of file + nil) + +(rt:deftest expintegral-e.4d + (let* ((x .5d0) + (e (exp-integral-e -2 x)) + (true (/ (* (exp (- x)) (+ (* x x x) (* 2 x x) (* 2 x))) + (expt x 4)))) + (check-accuracy 53 e true)) + nil) + +(rt:deftest expintegral-e.4q + (let* ((x #q.5) + (e (exp-integral-e -2 x)) + (true (/ (* (exp (- x)) (+ (* x x x) (* 2 x x) (* 2 x))) + (expt x 4)))) + (check-accuracy 210.8 e true)) + nil) ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 25 ++++++++++++++----------- rt-tests.lisp | 30 ++++++++++++++++++++++++++++++ 2 files changed, 44 insertions(+), 11 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 4b332ed2140e56c6fcaa689e6b8a48be36c988f6
by Raymond Toy 22 Mar '12

22 Mar '12

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 4b332ed2140e56c6fcaa689e6b8a48be36c988f6 (commit) from 9d3daf46c3e396941c8eb43209a45c0105217840 (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 4b332ed2140e56c6fcaa689e6b8a48be36c988f6 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Wed Mar 21 21:53:37 2012 -0700 Fix bug in psi for -n/2 for n odd which was causing an overflow. Add tests too. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index 02b5547..d49acf8 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -707,6 +707,8 @@ (realpart (ci z)) (ci z)))) +;; Array of values of the Bernoulli numbers. We only have enough for +;; the evaluation of the psi function. (defconstant bern-values (make-array 55 :initial-contents @@ -793,8 +795,15 @@ ;; formula to increase the argument and then apply the asymptotic formula. (cond ((minusp (realpart z)) - (- (psi (- 1 z)) - (* +pi+ (/ (tan (* +pi+ z)))))) + (let ((p (float +pi+ (realpart z)))) + (flet ((cot-pi (z) + ;; cot(%pi*z), car + (handler-case + (/ (tan (* p z))) + (division-by-zero () + (* 0 z))))) + (- (psi (- 1 z)) + (* p (cot-pi z)))))) (t (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon z) 10))))))) (m 0) diff --git a/rt-tests.lisp b/rt-tests.lisp index 9426c98..c5b6743 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1363,6 +1363,58 @@ (check-accuracy 212 s true)) nil) +(rt:deftest psi.1d + (let* ((z 1d0) + (p (psi z)) + (true (float (- +%gamma+) 1d0))) + (check-accuracy 52 p true)) + nil) + +(rt:deftest psi.1q + (let* ((z #q1) + (p (psi z)) + (true (- +%gamma+))) + (check-accuracy 208 p true)) + nil) + +(rt:deftest psi.2d + (let* ((z (float 4/3 1d0)) + (p (psi z)) + (true (- 3 + +%gamma+ + (/ +pi+ (* 2 (sqrt #q3))) + (* 1.5 (log #q3))))) + (check-accuracy 49.8 p true)) + nil) + +(rt:deftest psi.2d + (let* ((z (float 4/3 #q1)) + (p (psi z)) + (true (- 3 + +%gamma+ + (/ +pi+ (* 2 (sqrt #q3))) + (* 1.5 (log #q3))))) + (check-accuracy 205 p true)) + nil) + +(rt:deftest psi.3d + (let* ((z (float -1/2 1d0)) + (p (psi z)) + (true (- 2 + +%gamma+ + (log #q4)))) + (check-accuracy 48 p true)) + nil) + +(rt:deftest psi.3q + (let* ((z (float -1/2 #q1)) + (p (psi z)) + (true (- 2 + +%gamma+ + (log #q4)))) + (check-accuracy 204.1 p true)) + nil) + (rt:deftest expintegral-e.1d (let* ((z 1d0) (e (exp-integral-e 0 z)) ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 13 +++++++++++-- rt-tests.lisp | 52 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 2 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 9d3daf46c3e396941c8eb43209a45c0105217840
by Raymond Toy 22 Mar '12

22 Mar '12

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 9d3daf46c3e396941c8eb43209a45c0105217840 (commit) from 0d5870201359817c679921a2d740fdd1697469b2 (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 9d3daf46c3e396941c8eb43209a45c0105217840 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Wed Mar 21 18:58:02 2012 -0700 Add a few tests for exp-integral-e. More needed. diff --git a/rt-tests.lisp b/rt-tests.lisp index 33b2444..9426c98 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1363,3 +1363,36 @@ (check-accuracy 212 s true)) nil) +(rt:deftest expintegral-e.1d + (let* ((z 1d0) + (e (exp-integral-e 0 z)) + (true (/ (exp (- z)) z))) + (check-accuracy 53 e true)) + nil) + +(rt:deftest expintegral-e.1q + (let* ((z #q1) + (e (exp-integral-e 0 z)) + (true (/ (exp (- z)) z))) + (check-accuracy 212 e true)) + nil) + +(rt:deftest expintegral-e.2d + (let* ((z 15d0) + (e (exp-integral-e 0 z)) + (true (/ (exp (- z)) z))) + (check-accuracy 53 e true)) + nil) + +(rt:deftest expintegral-e.2q + (let* ((z #q15) + (e (exp-integral-e 0 z)) + (true (/ (exp (- z)) z))) + (check-accuracy 212 e true)) + nil) + +(rt:deftest expintegral-e.3d + (let* ((e (exp-integral-e 2 1d0)) + (true 0.14849550677592204791835999d0)) + (check-accuracy 47.5 e true)) + nil) \ No newline at end of file ----------------------------------------------------------------------- Summary of changes: rt-tests.lisp | 33 +++++++++++++++++++++++++++++++++ 1 files changed, 33 insertions(+), 0 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 0d5870201359817c679921a2d740fdd1697469b2
by Raymond Toy 22 Mar '12

22 Mar '12

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 0d5870201359817c679921a2d740fdd1697469b2 (commit) from 405df618a38d3b8ddaae691f865bbf068e931ac5 (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 0d5870201359817c679921a2d740fdd1697469b2 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Wed Mar 21 18:45:44 2012 -0700 Implement psi and fix exp-integral-e for integral values of v. Needs some more work. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index 48f8345..02b5547 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -479,17 +479,44 @@ (* (- 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)) (-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))) + when (< (abs term) (* (abs sum) eps)) + 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: @@ -679,3 +706,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))))) + ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 145 ++++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 files changed, 139 insertions(+), 6 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 405df618a38d3b8ddaae691f865bbf068e931ac5
by Raymond Toy 21 Mar '12

21 Mar '12

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 405df618a38d3b8ddaae691f865bbf068e931ac5 (commit) from c388f81713d7b2a483000d3cee1af030ed2c1cac (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 405df618a38d3b8ddaae691f865bbf068e931ac5 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Wed Mar 21 09:35:22 2012 -0700 Add Euler's constant and fix dependencies. oct.asd: o qd-const2 depends on qd-const. qd-const.lisp: o Add qd value for Euler's constant qd-const2.lisp: o Add qd-real value for Euler's constant. diff --git a/oct.asd b/oct.asd index ed7a586..3f8d70c 100644 --- a/oct.asd +++ b/oct.asd @@ -52,7 +52,7 @@ :depends-on ("qd" "qd-const")) (:file "qd-class" :depends-on ("qd-fun")) - (:file "qd-const2" :depends-on ("qd-class")) + (:file "qd-const2" :depends-on ("qd-class" "qd-const")) (:file "qd-methods" :depends-on ("qd-class")) (:file "qd-reader" diff --git a/qd-const.lisp b/qd-const.lisp index 628a520..56840b7 100644 --- a/qd-const.lisp +++ b/qd-const.lisp @@ -3752,6 +3752,15 @@ (%make-qd-d q0 q1 q2 q3)))) "A table of cos(k/1024) for k = 1 to 256") +(defconstant +qd-%gamma+ + (multiple-value-bind (q0 q1 q2 q3) + (renorm-5 (scale-float (float 5199096506725913 1.0d0) -53) + (scale-float (float -6416270704492459 1.0d0) -110) + (scale-float (float -5430032177773649 1.0d0) -164) + (scale-float (float 7163319935831933 1.0d0) -218) + (scale-float (float -6507777405222912 1.0d0) -273)) + (%make-qd-d q0 q1 q2 q3)) + "Quad-double value of Euler's constant") #|| diff --git a/qd-const2.lisp b/qd-const2.lisp index 2ad911d..e5e65e6 100644 --- a/qd-const2.lisp +++ b/qd-const2.lisp @@ -82,3 +82,7 @@ (defconstant +qd-real-one+ (make-instance 'qd-real :value (make-qd-d 1d0)) "QD-REAL representation of 1") + +(defconstant +%gamma+ + (make-instance 'qd-real :value octi::+qd-%gamma+) + "Euler's constant") \ No newline at end of file ----------------------------------------------------------------------- Summary of changes: oct.asd | 2 +- qd-const.lisp | 9 +++++++++ qd-const2.lisp | 4 ++++ 3 files changed, 14 insertions(+), 1 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0