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 a44327a2532c5afc3cbcb1194aa1f15b6dccd2cb (commit) via 7fa2a4c33e7542f3538bfcfd15cba833120b0daf (commit) via bed5763d1250d8b3b3793861f6e17449c49dd15b (commit) via 70745f486e3ee012f02c080de0cf9daf00f37375 (commit) from cc23f1098858d652c66728ea85c1c331bfb0fbb8 (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 a44327a2532c5afc3cbcb1194aa1f15b6dccd2cb Author: Raymond Toy <toy.raymond@gmail.com> Date: Thu Mar 24 23:01:51 2011 -0400 Update required accuracy for elliptic-pi.n0.q. diff --git a/rt-tests.lisp b/rt-tests.lisp index ecb401a..c0b59fb 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1037,7 +1037,7 @@ for epi = (elliptic-pi n phi 0) for true = (/ (atan (* (tan phi) (sqrt (- 1 n)))) (sqrt (- 1 n))) - for result = (check-accuracy 208 epi true) + for result = (check-accuracy 204 epi true) unless (eq nil result) append (list (list (list k n phi) result))) nil) commit 7fa2a4c33e7542f3538bfcfd15cba833120b0daf Author: Raymond Toy <toy.raymond@gmail.com> Date: Thu Mar 24 23:01:08 2011 -0400 Add DOMAIN-ERROR condition; fix bug in FRESNEL-S-SERIES. qd-methods.lisp: o Define DOMAIN-ERROR condition to allow signaling errors for incorrect domains. qd-gamma.lisp: o Signal domain error in CF-INCOMPLETE-GAMMA-TAIL if necessary. o Fix bug in FRESNEL-S-SERIES. We were comparing a real against a complex. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index 74a9032..78b5a99 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -235,6 +235,11 @@ ;; ;; See http://functions.wolfram.com/06.06.10.0003.01 (defun cf-incomplete-gamma-tail (a z) + (when (and (zerop (imagpart z)) (minusp (realpart z))) + (error 'domain-error + :function-name 'cf-incomplete-gamma-tail + :format-arguments (list 'z z) + :format-control "Argument ~S should not be on the negative real axis: ~S")) (/ (* (expt z a) (exp (- z))) (let ((z-a (- z a))) @@ -391,7 +396,7 @@ (sum 0) (term pi/2)) (loop for k2 from 0 by 2 - until (< (abs term) (* eps sum)) + until (< (abs term) (* eps (abs sum))) do (incf sum (/ term (+ 3 k2 k2))) (setf term (/ (* term factor) diff --git a/qd-methods.lisp b/qd-methods.lisp index 29b296d..e4df131 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -1123,3 +1123,14 @@ the same precision as the argument. The argument can be complex.")) (defmethod float-pi ((z qd-complex)) +pi+) + +(define-condition domain-error (simple-error) + ((function-name :accessor condition-function-name + :initarg :function-name)) + (:report (lambda (condition stream) + (format stream "Domain Error for function ~S:~&" + (condition-function-name condition)) + (pprint-logical-block (stream nil :per-line-prefix " ") + (apply #'format stream + (simple-condition-format-control condition) + (simple-condition-format-arguments condition)))))) \ No newline at end of file commit bed5763d1250d8b3b3793861f6e17449c49dd15b Author: Raymond Toy <toy.raymond@gmail.com> Date: Thu Mar 24 14:48:55 2011 -0400 Fix incomplete-gamma-tail for negative reals, use series for Fresnel S for small arg and update tests. qd-gamma.lisp: o INCOMPLETE-GAMMA-TAIL was hanging for arguments on the negative real axis. Use INCOMPLETE-GAMMA in this case too. o Add the series expansion for Fresnel S and use it for evaluating it for small arguments. We were losing accuracy with the existing algorithm. rt-tests.lisp: o Update thresholds for elliptic-pi-n0.d, elliptic-pi.n2.q, theta3.1.d. o Fix typo in test name. gamma-incomplete-tail.1.q should have been 2.q. o Add tests for gamma-incomplete-tail for arguments on the negative real axis. o Add tests for Fresnel S. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index c0cacaa..74a9032 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -288,10 +288,17 @@ (let* ((prec (float-contagion a z)) (a (apply-contagion a prec)) (z (apply-contagion z prec))) - (if (and (realp a) (realp z)) + (if (and (zerop (imagpart a)) + (zerop (imagpart z))) ;; For real values, we split the result to compute either the - ;; tail directly or from gamma(a) - incomplete-gamma - (if (> (abs z) (abs (- a 1))) + ;; tail directly from the continued fraction or from gamma(a) + ;; - incomplete-gamma. The continued fraction doesn't + ;; converge on the negative real axis, so we can't use that + ;; there. And accuracy appears to be better if z is "small". + ;; We take this to mean |z| < |a-1|. Note that |a-1| is the + ;; peak of the integrand. + (if (and (> (abs z) (abs (- a 1))) + (not (minusp (realpart z)))) (cf-incomplete-gamma-tail a z) (- (gamma a) (incomplete-gamma a z))) (cf-incomplete-gamma-tail a z)))) @@ -364,27 +371,90 @@ (* (expt z (- v 1)) (incomplete-gamma-tail (- 1 v) z))) +;; Series for Fresnel S +;; +;; S(z) = z^3*sum((%pi/2)^(2*k+1)(-z^4)^k/(2*k+1)!/(4*k+3), k, 0, inf) +;; +;; Compute as +;; +;; S(z) = z^3*sum(a(k)/(4*k+3), k, 0, inf) +;; +;; where +;; +;; a(k+1) = -a(k) * (%pi/2)^2 * z^4 / (2*k+2) / (2*k+3) +;; +;; a(0) = %pi/2. +(defun fresnel-s-series (z) + (let* ((pi/2 (* 1/2 (float-pi z))) + (factor (- (* (expt z 4) pi/2 pi/2))) + (eps (epsilon z)) + (sum 0) + (term pi/2)) + (loop for k2 from 0 by 2 + until (< (abs term) (* eps sum)) + do + (incf sum (/ term (+ 3 k2 k2))) + (setf term (/ (* term factor) + (* (+ k2 2) + (+ k2 3))))) + (* sum (expt z 3)))) + (defun fresnel-s (z) "Fresnel S: S(z) = integrate(sin(%pi*t^2/2), t, 0, z) " - (let ((sqrt-pi (sqrt (float-pi z)))) + (let ((prec (float-contagion z)) + (sqrt-pi (sqrt (float-pi z)))) (flet ((fs (z) ;; Wolfram gives ;; ;; S(z) = (1+%i)/4*(erf(c*z) - %i*erf(conjugate(c)*z)) ;; ;; where c = sqrt(%pi)/2*(1+%i). - (* #c(1/4 1/4) - (- (erf (* #c(1/2 1/2) sqrt-pi z)) - (* #c(0 1) - (erf (* #c(1/2 -1/2) sqrt-pi z))))))) - (if (realp z) - ;; FresnelS is real for a real argument. And it is odd. - (if (minusp z) - (- (realpart (fs (- z)))) - (realpart (fs z))) - (fs z))))) + ;; + ;; But for large z, we should use erfc. Then + ;; S(z) = 1/2 - (1+%i)/4*(erfc(c*z) - %i*erfc(conjugate(c)*z)) + (if (and t (> (abs z) 2)) + (- 1/2 + (* #c(1/4 1/4) + (- (erfc (* #c(1/2 1/2) sqrt-pi z)) + (* #c(0 1) + (erfc (* #c(1/2 -1/2) sqrt-pi z)))))) + (* #c(1/4 1/4) + (- (erf (* #c(1/2 1/2) sqrt-pi z)) + (* #c(0 1) + (erf (* #c(1/2 -1/2) sqrt-pi z))))))) + (rfs (z) + ;; When z is real, recall that erf(conjugate(z)) = + ;; conjugate(erf(z)). Then + ;; + ;; S(z) = 1/2*(realpart(erf(c*z)) - imagpart(erf(c*z))) + ;; + ;; But for large z, we should use erfc. Then + ;; + ;; S(z) = 1/2 - 1/2*(realpart(erfc(c*z)) - imagpart(erf(c*z))) + (if (> (abs z) 2) + (let ((s (erfc (* #c(1/2 1/2) sqrt-pi z)))) + (- 1/2 + (* 1/2 (- (realpart s) (imagpart s))))) + (let ((s (erf (* #c(1/2 1/2) sqrt-pi z)))) + (* 1/2 (- (realpart s) (imagpart s))))))) + ;; For small z, the erf terms above suffer from subtractive + ;; cancellation. So use the series in this case. Some simple + ;; tests were done to determine that for double-floats we want + ;; to use the series for z < 1 to give max accuracy. For + ;; qd-real, the above formula is good enough for z > 1d-5. + (if (< (abs z) (ecase prec + (single-float 1.5f0) + (double-float 1d0) + (qd-real #q1))) + (fresnel-s-series z) + (if (realp z) + ;; FresnelS is real for a real argument. And it is odd. + (if (minusp z) + (- (rfs (- z))) + (rfs z)) + (fs z)))))) (defun fresnel-c (z) "Fresnel C: diff --git a/rt-tests.lisp b/rt-tests.lisp index 5f4cdd5..ecb401a 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1000,7 +1000,7 @@ for epi = (elliptic-pi n phi 0) for true = (/ (atan (* (tan phi) (sqrt (- 1 n)))) (sqrt (- 1 n))) - for result = (check-accuracy 47.5 epi true) + for result = (check-accuracy 46.5 epi true) unless (eq nil result) append (list (list (list k n phi) result))) nil) @@ -1059,7 +1059,7 @@ for epi = (elliptic-pi n phi 0) for true = (/ (atanh (* (tan phi) (sqrt (- n 1)))) (sqrt (- n 1))) - for result = (check-accuracy 206 epi true) + for result = (check-accuracy 202 epi true) ;; Not sure if this formula holds when atanh gives a complex ;; result. Wolfram doesn't say when (and (not (complexp true)) result) @@ -1075,7 +1075,7 @@ for m = (random 1d0) for t3 = (elliptic-theta-3 0 (elliptic-nome m)) for true = (sqrt (/ (* 2 (elliptic-k m)) (float-pi m))) - for result = (check-accuracy 51 t3 true) + for result = (check-accuracy 50.5 t3 true) when result append (list (list (list k m) result))) nil) @@ -1213,15 +1213,74 @@ nil) (rt:deftest gamma-incomplete-tail.1.q - (let* ((z 5d0) + (let* ((z #q5) (gi (incomplete-gamma-tail 2 z)) (true (* (+ z 1) (exp (- z))))) - (check-accuracy 212 gi true)) + (check-accuracy 207 gi true)) nil) -(rt:deftest gamma-incomplete-tail.1.q +(rt:deftest gamma-incomplete-tail.2.q (let* ((z #q(1 5)) (gi (incomplete-gamma-tail 2 z)) (true (* (+ z 1) (exp (- z))))) (check-accuracy 206 gi true)) nil) + +(rt:deftest gamma-incomplete-tail.3.d + (let* ((z -5d0) + (gi (incomplete-gamma-tail 2 z)) + (true (* (+ z 1) (exp (- z))))) + (check-accuracy 52 gi true)) + nil) + +(rt:deftest gamma-incomplete-tail.3.q + (let* ((z #q-5) + (gi (incomplete-gamma-tail 2 z)) + (true (* (+ z 1) (exp (- z))))) + (check-accuracy 206 gi true)) + nil) + + +;; Fresnel integrals. +;; +;; For x small, Fresnel +;; +;; S(z) = %pi/6*z^3*(1 - %pi^2*z^4/56 + %pi^4*z^8/2040 - ...) +;; +(defun fresnel-s-series (z) + (let* ((fpi (float-pi z)) + (z^3 (expt z 3)) + (z^4 (* z^3 z))) + (* fpi 1/6 z^3 + (+ 1 (/ (* fpi fpi z^4) + -56) + (/ (* (expt fpi 4) (expt z^4 2)) + 7040))))) + +(rt:deftest fresnel-s.1d + (let* ((z 1d-3) + (s (fresnel-s z)) + (true (fresnel-s-series z))) + (check-accuracy 52 s true)) + nil) + +(rt:deftest fresnel-s.2d + (let* ((z #c(1d-3 1d-3)) + (s (fresnel-s z)) + (true (fresnel-s-series z))) + (check-accuracy 52 s true)) + nil) + +(rt:deftest fresnel-s.1q + (let* ((z #q1q-20) + (s (fresnel-s z)) + (true (fresnel-s-series z))) + (check-accuracy 212 s true)) + nil) + +(rt:deftest fresnel-s.2q + (let* ((z #q(1q-3 1q-3)) + (s (fresnel-s z)) + (true (fresnel-s-series z))) + (check-accuracy 212 s true)) + nil) commit 70745f486e3ee012f02c080de0cf9daf00f37375 Author: Raymond Toy <toy.raymond@gmail.com> Date: Thu Mar 24 14:25:53 2011 -0400 qd-gamma depends on qd-reader. diff --git a/oct.asd b/oct.asd index e4011a9..7df866c 100644 --- a/oct.asd +++ b/oct.asd @@ -66,7 +66,7 @@ (:file "qd-theta" :depends-on ("qd-methods" "qd-reader")) (:file "qd-gamma" - :depends-on ("qd-methods")) + :depends-on ("qd-methods" "qd-reader")) )) (defmethod perform ((op test-op) (c (eql (find-system :oct)))) ----------------------------------------------------------------------- Summary of changes: oct.asd | 2 +- qd-gamma.lisp | 103 +++++++++++++++++++++++++++++++++++++++++++++++------- qd-methods.lisp | 11 ++++++ rt-tests.lisp | 73 +++++++++++++++++++++++++++++++++++---- 4 files changed, 167 insertions(+), 22 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats