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 94d15953c5116999ac15377515b6c7e5ae430ea7 (commit) from de65a54b545af76b6df16f3de36e54763ce97235 (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 94d15953c5116999ac15377515b6c7e5ae430ea7 Author: Raymond Toy <toy.raymond@gmail.com> Date: Thu Apr 25 20:11:27 2013 -0700 Add continued fractions for erf, erfc, w and Dawson's integral. qd-gamma.lisp:: Add continued fractions rt-tests.lisp:: Add a test for erfc diff --git a/qd-gamma.lisp b/qd-gamma.lisp index 40a45eb..11e9822 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -275,6 +275,107 @@ #'(lambda (n) (* 4 n z2))))))) +;; From Cuyt +;; +;; sqrt(%pi)*z*exp(z^2)*erf(z) = K +;; +;; where K is the continued fraction with terms F[m]*z^2/(1 + G[m]*z^2) +;; with F[1] = 2 and F[m] = 4*(m-1)/(2*m-3)/(2*m-1) and G[m] = -2/(2*m-1) +;; +(defun cf-erf (z) + (let ((z2 (* z z))) + (* (/ (exp (- z2)) + (sqrt (float-pi z)) + z) + (lentz #'(lambda (n) + (if (zerop n) + (float 0 (realpart z)) + (1+ (/ (* -2 z2) + (+ n n -1))))) + #'(lambda (n) + (if (= n 1) + (* 2 z2) + (/ (* z2 4 (- n 1)) + (* (+ n n -3) (+ n n -1))))))))) + +;; From the above, we also have Dawson's integral: +;; +;; exp(-z^-2)*integrate(exp(t^2), t, 0, z) = %i*sqrt(%pi)/2*exp(-z^2)*erf(-%i*z); +;; +;; -2*z*exp(-z^2)*integrate(exp(t^2), t, 0, z) = K +;; +;; with K = -F[m)*z^2/(1 - G[m]*z^2), where F[m] and G[m] are as above. +;; +;; Also erf(-%i*z) = dawson(z) * 2*exp(-z^2)/(*%i*sqrt(%pi)) +(defun cf-dawson (z) + (let ((z2 (* z z))) + (/ (lentz #'(lambda (n) + (if (zerop n) + (float 0 (realpart z)) + (- 1 (/ (* -2 z2) + (+ n n -1))))) + #'(lambda (n) + (if (= n 1) + (* -2 z2) + (/ (* z2 -4 (- n 1)) + (* (+ n n -3) (+ n n -1)))))) + (* -2 z)))) + +;; erfc(z) = z/sqrt(%pi)*exp(-z^2)*K +;; +;; where K is the continued fraction with a[1] = 1, a[m] = (m-1)/2, +;; for m >= 2 and b[0] = 0, b[2*m+1] = z^2, b[2*m] = 1. +;; +;; This is valid only if Re(z) > 0. + +(defun cf-erfc (z) + (let ((z2 (* z z)) + (zero (float 0 (realpart z))) + (one (float 1 (realpart z)))) + (* (exp (- z2)) + z + (/ (sqrt (float-pi (realpart z)))) + (lentz #'(lambda (n) + (if (zerop n) + zero + (if (evenp n) + one + z2))) + #'(lambda (n) + (if (= n 1) + one + (/ (- n 1) 2))))))) + +;; w(z) = exp(-z^2)*erfc(-%i*z) +;; +;; = -%i*z/sqrt(%pi)*K +;; +;; where K is the continued fraction with a[n] the same as for erfc +;; and b[0] = 0, b[2*m+1] = -z^2, b[2*m] = 1. +;; +;; This is valid only if Im(z) > 0. We can use the following +;; identities: +;; +;; w(-z) = 2*exp(-z^2) - w(z) +;; w(conj(z)) = conj(w(-z)) + +(defun cf-w (z) + (let ((z2 (* z z)) + (zero (float 0 (realpart z))) + (one (float 1 (realpart z)))) + (* #c(0 -1) + z + (/ (sqrt (float-pi (realpart z)))) + (lentz #'(lambda (n) + (if (zerop n) + zero + (if (evenp n) + one + (- z2)))) + #'(lambda (n) + (if (= n 1) + one + (/ (- n 1) 2))))))) ;; Tail of the incomplete gamma function: ;; integrate(x^(a-1)*exp(-x), x, z, inf) ;; diff --git a/rt-tests.lisp b/rt-tests.lisp index af1d00c..caf568f 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1729,3 +1729,8 @@ when result append (list (list (list k x) result))) nil) + +(rt:deftest erfc + (check-accuracy 210 (erfc #q-4) + #q1.9999999845827420997199811478403265131159514278547464108088316570950057869589732) + nil) ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 101 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rt-tests.lisp | 5 +++ 2 files changed, 106 insertions(+) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats