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 c6f0a5ef00138c7f504432e0d0f485c23328274a (commit) via 80be2f3a917f7f429012985cd3a003d77de22ce7 (commit) from d2650578689f51edbdc8f0c588fcf330134490c6 (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 c6f0a5ef00138c7f504432e0d0f485c23328274a Author: Raymond Toy <toy.raymond@gmail.com> Date: Mon Mar 14 20:47:19 2011 -0400 Add docstrings; simplify elliptic-c, elliptic-ec. diff --git a/qd-elliptic.lisp b/qd-elliptic.lisp index 01a9e21..385f00f 100644 --- a/qd-elliptic.lisp +++ b/qd-elliptic.lisp @@ -413,7 +413,13 @@ ;; Elliptic integral of the first kind. This is computed using ;; Carlson's Rf function: ;; -;; F(phi, m) = sin(phi) * Rf(cos(phi)^2, 1 - m*sin(phi)^2, 1) +;; F(phi, m) = sin(phi) * Rf(cos(phi)^2, 1 - m*sin(phi)^2, 1) +;; +;; Also, DLMF 19.25.5 says +;; +;; F(phi, k) = Rf(c-1, c-k^2, c) +;; +;; where c = csc(phi)^2. (defun elliptic-f (x m) "Incomplete Elliptic integral of the first kind: @@ -485,6 +491,30 @@ ;; ;; E(phi, m) = integrate(sqrt(1-m*sin(x)^2), x, 0, phi) ;; +;; Carlson says +;; +;; E(phi, k) = sin(phi) * Rf(cos(phi)^2, 1-k^2*sin(phi)^2, 1) +;; - (k^2/3)*sin(phi)^3 * Rd(cos(phi)^2, 1 - k^2*sin(phi)^2, 1) +;; +;; But note that DLMF 19.25.9 says +;; +;; E(phi, k) = Rf(c-1, c - k^2, c) - k^2/3*Rd(c-1, c-k^2, c) +;; +;; where c = csc(phi)^2. Also DLMF 19.25.10 has +;; +;; E(phi, k) = k1^2*Rf(c-1,c-k^2,c) + k^/3*k1^2*Rd(c-1, c, c-k^2) +;; + k^2*sqrt((c-1)/(c*(c-k^2))) +;; +;; where k1^2 = 1-k^2 and c > k^2. Finally, DLMF 19.25.11 says +;; +;; E(phi, k) = -k1^2/3*Rd(c-k^2,c,c-1) + sqrt((c-k^2)/(c*(c-1))) +;; +;; for phi /= pi/2. +;; +;; One possible advantage is that all terms on the rhs are positive +;; for 19.25.9 if k^2 <= 0; 19.25.10 if 0 <= k^2 <= 1; 19.25.11 if 1 +;; <= k^2 <= c. It might be beneficial to use this idea so that we +;; don't have subtractive cancellation. (defun elliptic-e (phi m) "Incomplete elliptic integral of the second kind: @@ -499,21 +529,26 @@ E(phi, m) = integrate(sqrt(1-m*sin(x)^2), x, 0, phi)" ;; A&S 17.4.25 (sin phi)) (t - (let* ((sin-phi (sin phi)) - (cos-phi (cos phi)) - (k (sqrt m)) - (y (* (- 1 (* k sin-phi)) - (+ 1 (* k sin-phi))))) - (- (* sin-phi - (carlson-rf (* cos-phi cos-phi) y 1)) - (* (/ m 3) - (expt sin-phi 3) - (carlson-rd (* cos-phi cos-phi) y 1)))))))) + ;; For quad-doubles, it's significantly faster to compute + ;; cis(phi) than to compute sin and cos separately. + (multiple-value-bind (cos-phi sin-phi) + (let ((cis (cis phi))) + (values (realpart cis) (imagpart cis))) + (let ((y (- 1 (* m sin-phi sin-phi)))) + (- (* sin-phi + (carlson-rf (* cos-phi cos-phi) y 1)) + (* (/ m 3) + (expt sin-phi 3) + (carlson-rd (* cos-phi cos-phi) y 1))))))))) ;; Complete elliptic integral of second kind. ;; ;; E(phi) = integrate(sqrt(1-m*sin(x)^2), x, 0, %pi/2) ;; +;; Carlson says +;; +;; E(k) = Rf(0, 1-k^2,1) - (k^2/3)*Rd(0,1-k^2,1) +;; (defun elliptic-ec (m) "Complete elliptic integral of the second kind: @@ -525,12 +560,10 @@ E(m) = integrate(sqrt(1-m*sin(x)^2), x, 0, %pi/2)" ;; A&S 17.4.25 (float 1 m)) (t - (let* ((k (sqrt m)) - (y (* (- 1 k) - (+ 1 k)))) - (- (carlson-rf 0 y 1) + (let* ((m1 (- 1 m))) + (- (carlson-rf 0 m1 1) (* (/ m 3) - (carlson-rd 0 y 1))))))) + (carlson-rd 0 m1 1))))))) ;; Carlson's Rc function. ;; commit 80be2f3a917f7f429012985cd3a003d77de22ce7 Author: Raymond Toy <toy.raymond@gmail.com> Date: Mon Mar 14 20:45:28 2011 -0400 Export functions, add docstrings. qd-package.lisp: o Export the elliptic theta functions and the Carlson elliptic integrals. qd-theta.lisp: o Add simple docstrings for elliptic theta functions. o Add ELLIPTIC-THETA function. diff --git a/qd-package.lisp b/qd-package.lisp index 9587038..b102435 100644 --- a/qd-package.lisp +++ b/qd-package.lisp @@ -286,7 +286,15 @@ #:elliptic-k #:elliptic-f #:elliptic-e - #:elliptic-ec) + #:elliptic-ec + #:carlson-rd + #:carlson-rf + #:carlson-rj + #:elliptic-theta-1 + #:elliptic-theta-2 + #:elliptic-theta-3 + #:elliptic-theta-4 + #:elliptic-theta) ;; Constants (:export #:+pi+ #:+pi/2+ diff --git a/qd-theta.lisp b/qd-theta.lisp index caae788..fb41b82 100644 --- a/qd-theta.lisp +++ b/qd-theta.lisp @@ -70,6 +70,9 @@ ;; [ 0 0 1 ] (defun elliptic-theta-1 (z q) + "Elliptic theta function 1 + + theta1(z, q) = 2*q^(1/4)*sum((-1)^n*q^(n*(n+1))*sin((2*n+1)*z), n, 0, inf)" (let* ((precision (float-contagion z q)) (z (apply-contagion z precision)) (q (apply-contagion q precision)) @@ -97,6 +100,9 @@ ;; k = 1 [ ] ;; [ 0 0 1 ] (defun elliptic-theta-3 (z q) + "Elliptic theta function 3 + + theta3(z, q) = 1 + 2 * sum(q^(n^2)*cos(2*n*z), n, 1, inf)" (let* ((precision (float-contagion z q)) (z (apply-contagion z precision)) (q (apply-contagion q precision)) @@ -115,6 +121,9 @@ ;; theta[2](z,q) = theta[1](z+%pi/2, q) (defun elliptic-theta-2 (z q) + "Elliptic theta function 2 + + theta2(z, q) = 2*q^(1/4)*sum(q^(n*(n+1))*cos((2*n+1)*z), n, 0, inf)" (let* ((precision (float-contagion z q)) (z (apply-contagion z precision)) (q (apply-contagion q precision))) @@ -122,14 +131,26 @@ ;; theta[4](z,q) = theta[3](z+%pi/2,q) (defun elliptic-theta-4 (z q) + "Elliptic theta function 4 + + theta4(z, q) = 1 + 2*sum((-1)^n*q^(n^2)*cos(2*n*z), n, 1, inf)" (let* ((precision (float-contagion z q)) (z (apply-contagion z precision)) (q (apply-contagion q precision))) (elliptic-theta-3 (+ z (/ (float-pi z) 2)) q))) +(defun elliptic-theta (n z q) + "Elliptic Theta function n where n = 1, 2, 3, or 4." + (ecase n + (1 (elliptic-theta-1 z q)) + (2 (elliptic-theta-2 z q)) + (3 (elliptic-theta-3 z q)) + (4 (elliptic-theta-4 z q)))) + ;; The nome, q, is given by q = exp(-%pi*K'/K) where K and %i*K' are ;; the quarter periods. (defun elliptic-nome (m) + "Compute the elliptic nome, q, from the parameter m" (exp (- (/ (* (float-pi m) (elliptic-k (- 1 m))) (elliptic-k m))))) ----------------------------------------------------------------------- Summary of changes: qd-elliptic.lisp | 65 ++++++++++++++++++++++++++++++++++++++++------------- qd-package.lisp | 10 +++++++- qd-theta.lisp | 21 +++++++++++++++++ 3 files changed, 79 insertions(+), 17 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats