[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. a44327a2532c5afc3cbcb1194aa1f15b6dccd2cb
Raymond Toy
rtoy at common-lisp.net
Fri Mar 25 03:02:05 UTC 2011
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 at 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 at 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 at 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 at 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
More information about the oct-cvs
mailing list