[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. c6f0a5ef00138c7f504432e0d0f485c23328274a
Raymond Toy
rtoy at common-lisp.net
Tue Mar 15 00:52:30 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 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 at 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 at 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
More information about the oct-cvs
mailing list