[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 5bebc8d76a83b85cdba6d803d026aae741b6e7c8
Raymond Toy
rtoy at common-lisp.net
Fri Mar 11 04:22:53 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 5bebc8d76a83b85cdba6d803d026aae741b6e7c8 (commit)
via 5f6c628014268b683e61baae5470a56b078d1c16 (commit)
from c239a8514ff5f1fc96d33a4c581cbaec834f6641 (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 5bebc8d76a83b85cdba6d803d026aae741b6e7c8
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Thu Mar 10 23:05:08 2011 -0500
Test carlson-rf and carlson-rd.
diff --git a/rt-tests.lisp b/rt-tests.lisp
index 82ca0b7..85f9b95 100644
--- a/rt-tests.lisp
+++ b/rt-tests.lisp
@@ -833,3 +833,34 @@
(val (jacobi-dn ek #q.5)))
(check-accuracy 212 val true))
nil)
+
+(rt:deftest oct.carlson-rf.1d
+ ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
+ ;; = 1/4*beta(1/2,1/2)
+ ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
+ (let ((rf (carlson-rf 0d0 2d0 1d0))
+ (true 1.31102877714605990523241979494d0))
+ (check-accuracy 53 rf true))
+ nil)
+
+(rt:deftest oct.carlson-rf.1q
+ ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
+ (let ((rf (carlson-rf #q0 #q2 #q1))
+ (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0))
+ (check-accuracy 212 rf true))
+ nil)
+
+(rt:deftest oct.carlson-rd.1d
+ ;; Rd(0,2,1) = 3*integrate(s^2/sqrt(1-s^4), s, 0 ,1)
+ ;; = 3*beta(3/4,1/2)/4
+ ;; = 3*sqrt(%pi)*gamma(3/4)/gamma(1/4)
+ (let ((rd (carlson-rd 0d0 2d0 1d0))
+ (true 1.7972103521033883d0))
+ (check-accuracy 51 rd true))
+ nil)
+
+(rt:deftest oct.carlson-rd.1q
+ (let ((rd (carlson-rd #q0 #q2 #q1))
+ (true #q1.797210352103388311159883738420485817340818994823477337395512429419599q0))
+ (check-accuracy 212 rd true))
+ nil)
commit 5f6c628014268b683e61baae5470a56b078d1c16
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Thu Mar 10 23:04:28 2011 -0500
Add elliptic integrals of the second kind.
diff --git a/qd-elliptic.lisp b/qd-elliptic.lisp
index 63c3687..6c8070a 100644
--- a/qd-elliptic.lisp
+++ b/qd-elliptic.lisp
@@ -321,6 +321,73 @@
(* -3/44 ee2 ee3))))
(/ s (sqrt an)))))
+;; rd(x,y,z) = integrate(3/2*(t+x)^(-1/2)*(t+y)^(-1/2)*(t+z)^(-3/2), t, 0, inf)
+;;
+;; E(K) = rf(0, 1-K^2, 1) - (K^2/3)*rd(0,1-K^2,1)
+;;
+;; B = integrate(s^2/sqrt(1-s^4), s, 0 ,1)
+;; = beta(3/4,1/2)/4
+;; = sqrt(%pi)*gamma(3/4)/gamma(1/4)
+;; = 1/3*rd(0,2,1)
+(defun carlson-rd (x y z)
+ "Compute Carlson's Rd function:
+
+ Rd(x,y,z) = integrate(3/2*(t+x)^(-1/2)*(t+y)^(-1/2)*(t+z)^(-3/2), t, 0, inf)"
+ (let* ((xn x)
+ (yn y)
+ (zn z)
+ (a (/ (+ xn yn (* 3 zn)) 5))
+ (epslon (/ (max (abs (- a xn))
+ (abs (- a yn))
+ (abs (- a zn)))
+ (errtol x y z)))
+ (an a)
+ (sigma 0)
+ (power4 1)
+ (n 0)
+ xnroot ynroot znroot lam)
+ (loop while (> (* power4 epslon) (abs an))
+ do
+ (setf xnroot (sqrt xn))
+ (setf ynroot (sqrt yn))
+ (setf znroot (sqrt zn))
+ (setf lam (+ (* xnroot ynroot)
+ (* xnroot znroot)
+ (* ynroot znroot)))
+ (setf sigma (+ sigma (/ power4
+ (* znroot (+ zn lam)))))
+ (setf power4 (* power4 1/4))
+ (setf xn (* (+ xn lam) 1/4))
+ (setf yn (* (+ yn lam) 1/4))
+ (setf zn (* (+ zn lam) 1/4))
+ (setf an (* (+ an lam) 1/4))
+ (incf n))
+ ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26
+ (let* ((xndev (/ (* (- a x) power4) an))
+ (yndev (/ (* (- a y) power4) an))
+ (zndev (- (* (+ xndev yndev) 1/3)))
+ (ee2 (- (* xndev yndev) (* 6 zndev zndev)))
+ (ee3 (* (- (* 3 xndev yndev)
+ (* 8 zndev zndev))
+ zndev))
+ (ee4 (* 3 (- (* xndev yndev) (* zndev zndev)) zndev zndev))
+ (ee5 (* xndev yndev zndev zndev zndev))
+ (s (+ 1
+ (* -3/14 ee2)
+ (* 1/6 ee3)
+ (* 9/88 ee2 ee2)
+ (* -3/22 ee4)
+ (* -9/52 ee2 ee3)
+ (* 3/26 ee5)
+ (* -1/16 ee2 ee2 ee2)
+ (* 3/10 ee3 ee3)
+ (* 3/20 ee2 ee4)
+ (* 45/272 ee2 ee2 ee3)
+ (* -9/68 (+ (* ee2 ee5) (* ee3 ee4))))))
+ (+ (* 3 sigma)
+ (/ (* power4 s)
+ (expt an 3/2))))))
+
;; Complete elliptic integral of the first kind. This can be computed
;; from Carlson's Rf function:
;;
@@ -408,3 +475,55 @@
(* (- 1 (* k sin-x))
(+ 1 (* k sin-x)))
1)))))))
+
+;; Incomplete elliptic integral of the second kind.
+;;
+;; E(phi, m) = integrate(sqrt(1-m*sin(x)^2), x, 0, phi)
+;;
+(defun elliptic-e (phi m)
+ "Incomplete elliptic integral of the second kind:
+
+E(phi, m) = integrate(sqrt(1-m*sin(x)^2), x, 0, phi)"
+ (let* ((precision (float-contagion phi m))
+ (phi (coerce phi precision))
+ (m (coerce m precision)))
+ (cond ((= m 0)
+ ;; A&S 17.4.23
+ phi)
+ ((= m 1)
+ ;; 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 (coerce 1 precision)))
+ (* (/ m 3)
+ (expt sin-phi 3)
+ (carlson-rd (* cos-phi cos-phi) y (coerce 1 precision)))))))))
+
+;; Complete elliptic integral of second kind.
+;;
+;; E(phi) = integrate(sqrt(1-m*sin(x)^2), x, 0, %pi/2)
+;;
+(defun elliptic-ec (m)
+ "Complete elliptic integral of the second kind:
+
+E(m) = integrate(sqrt(1-m*sin(x)^2), x, 0, %pi/2)"
+ (let ((precision (float-contagion m)))
+ (cond ((= m 0)
+ ;; A&S 17.4.23
+ (/ (float-pi m) 2))
+ ((= m 1)
+ ;; A&S 17.4.25
+ (coerce 1 precision))
+ (t
+ (let* ((k (sqrt m))
+ (y (* (- 1 k)
+ (+ 1 k))))
+ (- (carlson-rf 0.0 y 1.0)
+ (* (/ m 3)
+ (carlson-rd 0.0 y 1.0))))))))
-----------------------------------------------------------------------
Summary of changes:
qd-elliptic.lisp | 119 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
rt-tests.lisp | 31 ++++++++++++++
2 files changed, 150 insertions(+), 0 deletions(-)
hooks/post-receive
--
OCT: A portable Lisp implementation for quad-double precision floats
More information about the oct-cvs
mailing list