Changeset 390a7483f5658fe802d5d239070cdaa573adf4a5

Timestamp:

03/12/11 10:35:18 (2 years ago)

Author:

Raymond Toy <toy.raymond@…>

Children:

5bd5df9360268fae90fc41fcdf86b728f8a54e86

Parents:

147fa2c7aa5e99988a7c3fb35800865d22efbc51

git-committer:

Raymond Toy <toy.raymond@…> (03/12/11 10:35:18)

Message:

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.

Files:

• qd-elliptic.lisp (modified) (2 diffs)
• rt-tests.lisp (modified) (1 diff)

qd-elliptic.lisp

@@ -409 +409 @@
          (/ (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
@@ -533 +532 @@
               (* (/ 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))))))))

rt-tests.lisp

@@ -930 +930 @@
       (check-accuracy 212 rf true))
   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)
+||#