[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 390a7483f5658fe802d5d239070cdaa573adf4a5
Raymond Toy
rtoy at common-lisp.net
Sat Mar 12 19:00:40 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 390a7483f5658fe802d5d239070cdaa573adf4a5 (commit)
from 147fa2c7aa5e99988a7c3fb35800865d22efbc51 (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 390a7483f5658fe802d5d239070cdaa573adf4a5
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Sat Mar 12 13:35:18 2011 -0500
Add prelimary support for integrals of the 3rd kind.
qd-elliptic.lisp:
o Clean up for unused variable in ELLIPTIC-K
o Add Carlson's Rj functions
o Implement elliptic-pi using Carlson's method.
rt-tests.lisp:
o Add many tests for elliptic-pi. Some tests pass, and some fail. The
failing tests are not enabled because I don't know if the failure is
because the test itself is wrong or if the integral is wrong.
diff --git a/qd-elliptic.lisp b/qd-elliptic.lisp
index 342290e..1b80f13 100644
--- a/qd-elliptic.lisp
+++ b/qd-elliptic.lisp
@@ -408,8 +408,7 @@
(cond ((= m 0)
(/ (float +pi+ m) 2))
(t
- (let ((precision (float-contagion m)))
- (carlson-rf 0 (- 1 m) 1)))))
+ (carlson-rf 0 (- 1 m) 1))))
;; Elliptic integral of the first kind. This is computed using
;; Carlson's Rf function:
@@ -532,3 +531,186 @@ E(m) = integrate(sqrt(1-m*sin(x)^2), x, 0, %pi/2)"
(- (carlson-rf 0 y 1)
(* (/ m 3)
(carlson-rd 0 y 1)))))))
+
+;; Carlson's Rc function.
+;;
+;; Some interesting identities:
+;;
+;; log(x) = (x-1)*rc(((1+x)/2)^2, x), x > 0
+;; asin(x) = x * rc(1-x^2, 1), |x|<= 1
+;; acos(x) = sqrt(1-x^2)*rc(x^2,1), 0 <= x <=1
+;; atan(x) = x * rc(1,1+x^2)
+;; asinh(x) = x * rc(1+x^2,1)
+;; acosh(x) = sqrt(x^2-1) * rc(x^2,1), x >= 1
+;; atanh(x) = x * rc(1,1-x^2), |x|<=1
+;;
+
+(defun carlson-rc (x y)
+ "Compute Carlson's Rc function:
+
+ Rc(x,y) = integrate(1/2*(t+x)^(-1/2)*(t+y)^(-1), t, 0, inf)"
+ (let* ((precision (float-contagion x y))
+ (yn (apply-contagion y precision))
+ (x (apply-contagion x precision))
+ xn z w a an pwr4 n epslon lambda sn s)
+ (cond ((and (zerop (imagpart yn))
+ (minusp (realpart yn)))
+ (setf xn (- x y))
+ (setf yn (- yn))
+ (setf z yn)
+ (setf w (sqrt (/ x xn))))
+ (t
+ (setf xn x)
+ (setf z yn)
+ (setf w 1)))
+ (setf a (/ (+ xn yn yn) 3))
+ (setf epslon (/ (abs (- a xn)) (errtol x y)))
+ (setf an a)
+ (setf pwr4 1)
+ (setf n 0)
+ (loop while (> (* epslon pwr4) (abs an))
+ do
+ (setf pwr4 (/ pwr4 4))
+ (setf lambda (+ (* 2 (sqrt xn) (sqrt yn)) yn))
+ (setf an (/ (+ an lambda) 4))
+ (setf xn (/ (+ xn lambda) 4))
+ (setf yn (/ (+ yn lambda) 4))
+ (incf n))
+ ;; c2=3/10,c3=1/7,c4=3/8,c5=9/22,c6=159/208,c7=9/8
+ (setf sn (/ (* pwr4 (- z a)) an))
+ (setf s (* sn sn (+ 3/10
+ (* sn (+ 1/7
+ (* sn (+ 3/8
+ (* sn (+ 9/22
+ (* sn (+ 159/208
+ (* sn 9/8))))))))))))
+ (/ (* w (+ 1 s))
+ (sqrt an))))
+
+(defun carlson-rj1 (x y z p)
+ (let* ((xn x)
+ (yn y)
+ (zn z)
+ (pn p)
+ (en (* (- pn xn)
+ (- pn yn)
+ (- pn zn)))
+ (sigma 0)
+ (power4 1)
+ (k 0)
+ (a (/ (+ xn yn zn pn pn) 5))
+ (epslon (/ (max (abs (- a xn))
+ (abs (- a yn))
+ (abs (- a zn))
+ (abs (- a pn)))
+ (errtol x y z p)))
+ (an a)
+ xnroot ynroot znroot pnroot lam dn)
+ (loop while (> (* power4 epslon) (abs an))
+ do
+ (setf xnroot (sqrt xn))
+ (setf ynroot (sqrt yn))
+ (setf znroot (sqrt zn))
+ (setf pnroot (sqrt pn))
+ (setf lam (+ (* xnroot ynroot)
+ (* xnroot znroot)
+ (* ynroot znroot)))
+ (setf dn (* (+ pnroot xnroot)
+ (+ pnroot ynroot)
+ (+ pnroot znroot)))
+ (setf sigma (+ sigma
+ (/ (* power4
+ (carlson-rc 1 (+ 1 (/ en (* dn dn)))))
+ dn)))
+ (setf power4 (* power4 1/4))
+ (setf en (/ en 64))
+ (setf xn (* (+ xn lam) 1/4))
+ (setf yn (* (+ yn lam) 1/4))
+ (setf zn (* (+ zn lam) 1/4))
+ (setf pn (* (+ pn lam) 1/4))
+ (setf an (* (+ an lam) 1/4))
+ (incf k))
+ (let* ((xndev (/ (* (- a x) power4) an))
+ (yndev (/ (* (- a y) power4) an))
+ (zndev (/ (* (- a z) power4) an))
+ (pndev (* -0.5 (+ xndev yndev zndev)))
+ (ee2 (+ (* xndev yndev)
+ (* xndev zndev)
+ (* yndev zndev)
+ (* -3 pndev pndev)))
+ (ee3 (+ (* xndev yndev zndev)
+ (* 2 ee2 pndev)
+ (* 4 pndev pndev pndev)))
+ (ee4 (* (+ (* 2 xndev yndev zndev)
+ (* ee2 pndev)
+ (* 3 pndev pndev pndev))
+ pndev))
+ (ee5 (* xndev yndev zndev pndev pndev))
+ (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))))))
+ (+ (* 6 sigma)
+ (/ (* power4 s)
+ (sqrt (* an an an)))))))
+
+(defun carlson-rj (x y z p)
+ "Compute Carlson's Rj function:
+
+ Rj(x,y,z,p) = integrate(3/2*(t+x)^(-1/2)*(t+y)^(-1/2)*(t+z)^(-1/2)*(t+p)^(-1), t, 0, inf)"
+ (let* ((precision (float-contagion x y z p))
+ (xn (apply-contagion x precision))
+ (yn (apply-contagion y precision))
+ (zn (apply-contagion z precision))
+ (p (apply-contagion p precision))
+ (qn (- p)))
+ (cond ((and (and (zerop (imagpart xn)) (>= (realpart xn) 0))
+ (and (zerop (imagpart yn)) (>= (realpart yn) 0))
+ (and (zerop (imagpart zn)) (>= (realpart zn) 0))
+ (and (zerop (imagpart qn)) (> (realpart qn) 0)))
+ (destructuring-bind (xn yn zn)
+ (sort (list xn yn zn) #'<)
+ (let* ((pn (+ yn (* (- zn yn) (/ (- yn xn) (+ yn qn)))))
+ (s (- (* (- pn yn) (carlson-rj1 xn yn zn pn))
+ (* 3 (carlson-rf xn yn zn)))))
+ (setf s (+ s (* 3 (sqrt (/ (* xn yn zn)
+ (+ (* xn zn) (* pn qn))))
+ (carlson-rc (+ (* xn zn) (* pn qn)) (* pn qn)))))
+ (/ s (+ yn qn)))))
+ (t
+ (carlson-rj1 x y z p)))))
+
+;; Elliptic integral of the third kind:
+;;
+;; (A&S 17.2.14)
+;;
+;; PI(n; phi|m) = integrate(1/sqrt(1-m*sin(x)^2)/(1-n*sin(x)^2), x, 0, phi)
+;;
+(defun elliptic-pi (n phi m)
+ "Compute elliptic integral of the third kind:
+
+ PI(n; phi|m) = integrate(1/sqrt(1-m*sin(x)^2)/(1-n*sin(x)^2), x, 0, phi)"
+ ;; Note: Carlson's DRJ has n defined as the negative of the n given
+ ;; in A&S.
+ (let* ((precision (float-contagion n phi m))
+ (n (apply-contagion n precision))
+ (phi (apply-contagion phi precision))
+ (m (apply-contagion m precision))
+ (nn (- n))
+ (sin-phi (sin phi))
+ (cos-phi (cos phi))
+ (k (sqrt m))
+ (k2sin (* (- 1 (* k sin-phi))
+ (+ 1 (* k sin-phi)))))
+ (- (* sin-phi (carlson-rf (expt cos-phi 2) k2sin 1))
+ (* (/ nn 3) (expt sin-phi 3)
+ (carlson-rj (expt cos-phi 2) k2sin 1
+ (+ 1 (* nn (expt sin-phi 2))))))))
diff --git a/rt-tests.lisp b/rt-tests.lisp
index 7635611..80d5c57 100644
--- a/rt-tests.lisp
+++ b/rt-tests.lisp
@@ -928,4 +928,143 @@
(let ((rf (carlson-rf 0 2 #q1q0))
(true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0))
(check-accuracy 212 rf true))
- nil)
\ No newline at end of file
+ nil)
+
+;; Elliptic integral of the third kind
+
+;; elliptic-pi(0,phi,m) = elliptic-f(phi, m)
+(rt:deftest oct.elliptic-pi.1d
+ (loop for k from 0 to 100
+ for phi = (random (/ pi 2))
+ for m = (random 1d0)
+ for epi = (elliptic-pi 0 phi m)
+ for ef = (elliptic-f phi m)
+ for result = (check-accuracy 53 epi ef)
+ unless (eq nil result)
+ append (list (list phi m) result))
+ nil)
+
+(rt:deftest oct.elliptic-pi.1q
+ (loop for k from 0 below 100
+ for phi = (random (/ +pi+ 2))
+ for m = (random #q1)
+ for epi = (elliptic-pi 0 phi m)
+ for ef = (elliptic-f phi m)
+ for result = (check-accuracy 53 epi ef)
+ unless (eq nil result)
+ append (list (list phi m) result))
+ nil)
+
+;; DLMF 19.6.3
+;;
+;; PI(n; pi/2 | 0) = pi/(2*sqrt(1-n))
+(rt:deftest oct.elliptic-pi.19.6.3.d
+ (loop for k from 0 below 100
+ for n = (random 1d0)
+ for epi = (elliptic-pi n (/ pi 2) 0)
+ for true = (/ pi (* 2 (sqrt (- 1 n))))
+ for result = (check-accuracy 49 epi true)
+ unless (eq nil result)
+ append (list (list (list k n) result)))
+ nil)
+
+(rt:deftest oct.elliptic-pi.19.6.3.q
+ (loop for k from 0 below 100
+ for n = (random #q1)
+ for epi = (elliptic-pi n (/ (float-pi n) 2) 0)
+ for true = (/ (float-pi n) (* 2 (sqrt (- 1 n))))
+ for result = (check-accuracy 210 epi true)
+ unless (eq nil result)
+ append (list (list (list k n) result)))
+ nil)
+
+#+nil
+(rt:deftest oct.elliptic-pi.19.6.2.d
+ (loop for k from 0 below 100
+ for n = (random 1d0)
+ for epi = (elliptic-pi (- n) (/ (float-pi n) 2) n)
+ for true = (+ (/ (float-pi n) 4 (sqrt (+ 1 (sqrt n))))
+ (/ (elliptic-k n) 2))
+ for result = (check-accuracy 53 epi true)
+ when result
+ append (list (list (list k n) result)))
+ nil)
+
+
+#||
+;; elliptic-pi(n, phi, 0) =
+;; atanh(sqrt(1-n)*tan(phi))/sqrt(1-n) n < 1
+;; atanh(sqrt(n-1)*tan(phi))/sqrt(n-1) n > 1
+;; tan(phi) n = 1
+(rt:deftest oct.elliptic-pi.n0.d
+ (loop for k from 0 below 100
+ for phi = (random (/ pi 2))
+ for n = (random 1d0)
+ for epi = (elliptic-pi n phi 0)
+ for true = (/ (atanh (* (tan phi) (sqrt (- 1 n))))
+ (sqrt (- 1 n)))
+ for result = (check-accuracy 53 epi true)
+ unless (eq nil result)
+ append (list (list (list k n phi) result)))
+ nil)
+
+(rt:deftest oct.elliptic-pi.n1.d
+ (loop for k from 0 below 100
+ for phi = (random (/ pi 2))
+ for epi = (elliptic-pi 0 phi 0)
+ for true = (tan phi)
+ for result = (check-accuracy 53 epi true)
+ unless (eq nil result)
+ append (list (list (list k phi) result)))
+ nil)
+
+(rt:deftest oct.elliptic-pi.n2.d
+ (loop for k from 0 below 100
+ for phi = (random (/ pi 2))
+ for n = (+ 1d0 (random 100d0))
+ for epi = (elliptic-pi n phi 0)
+ for true = (/ (atanh (* (tan phi) (sqrt (- n 1))))
+ (sqrt (- n 1)))
+ for result = (check-accuracy 52 epi true)
+ ;; Not sure if this formula holds when atanh gives a complex
+ ;; result. Wolfram doesn't say
+ when (and (not (complexp true)) result)
+ append (list (list (list k n phi) result)))
+ nil)
+
+(rt:deftest oct.elliptic-pi.n0.q
+ (loop for k from 0 below 100
+ for phi = (random (/ +pi+ 2))
+ for n = (random #q1)
+ for epi = (elliptic-pi n phi 0)
+ for true = (/ (atanh (* (tan phi) (sqrt (- 1 n))))
+ (sqrt (- 1 n)))
+ for result = (check-accuracy 212 epi true)
+ unless (eq nil result)
+ append (list (list (list k n phi) result)))
+ nil)
+
+(rt:deftest oct.elliptic-pi.n1.q
+ (loop for k from 0 below 100
+ for phi = (random (/ +pi+ 2))
+ for epi = (elliptic-pi 0 phi 0)
+ for true = (tan phi)
+ for result = (check-accuracy 212 epi true)
+ unless (eq nil result)
+ append (list (list (list k phi) result)))
+ nil)
+
+(rt:deftest oct.elliptic-pi.n2.q
+ (loop for k from 0 below 100
+ for phi = (random (/ +pi+ 2))
+ for n = (+ #q1 (random #q1))
+ for epi = (elliptic-pi n phi 0)
+ for true = (/ (atanh (* (tan phi) (sqrt (- n 1))))
+ (sqrt (- n 1)))
+ for result = (check-accuracy 209 epi true)
+ ;; Not sure if this formula holds when atanh gives a complex
+ ;; result. Wolfram doesn't say
+ when (and (not (complexp true)) result)
+ append (list (list (list k n phi) result)))
+ nil)
+||#
-----------------------------------------------------------------------
Summary of changes:
qd-elliptic.lisp | 186 +++++++++++++++++++++++++++++++++++++++++++++++++++++-
rt-tests.lisp | 141 ++++++++++++++++++++++++++++++++++++++++-
2 files changed, 324 insertions(+), 3 deletions(-)
hooks/post-receive
--
OCT: A portable Lisp implementation for quad-double precision floats
More information about the oct-cvs
mailing list