Commit 390a7483f5658fe802d5d239070cdaa573adf4a5

Authored by Raymond Toy
1 parent 147fa2c7

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.
qd-elliptic.lisp
... ... @@ -408,8 +408,7 @@
408 408 (cond ((= m 0)
409 409 (/ (float +pi+ m) 2))
410 410 (t
411   - (let ((precision (float-contagion m)))
412   - (carlson-rf 0 (- 1 m) 1)))))
  411 + (carlson-rf 0 (- 1 m) 1))))
413 412  
414 413 ;; Elliptic integral of the first kind. This is computed using
415 414 ;; Carlson's Rf function:
... ... @@ -532,3 +531,186 @@ E(m) = integrate(sqrt(1-m*sin(x)^2), x, 0, %pi/2)"
532 531 (- (carlson-rf 0 y 1)
533 532 (* (/ m 3)
534 533 (carlson-rd 0 y 1)))))))
  534 +
  535 +;; Carlson's Rc function.
  536 +;;
  537 +;; Some interesting identities:
  538 +;;
  539 +;; log(x) = (x-1)*rc(((1+x)/2)^2, x), x > 0
  540 +;; asin(x) = x * rc(1-x^2, 1), |x|<= 1
  541 +;; acos(x) = sqrt(1-x^2)*rc(x^2,1), 0 <= x <=1
  542 +;; atan(x) = x * rc(1,1+x^2)
  543 +;; asinh(x) = x * rc(1+x^2,1)
  544 +;; acosh(x) = sqrt(x^2-1) * rc(x^2,1), x >= 1
  545 +;; atanh(x) = x * rc(1,1-x^2), |x|<=1
  546 +;;
  547 +
  548 +(defun carlson-rc (x y)
  549 + "Compute Carlson's Rc function:
  550 +
  551 + Rc(x,y) = integrate(1/2*(t+x)^(-1/2)*(t+y)^(-1), t, 0, inf)"
  552 + (let* ((precision (float-contagion x y))
  553 + (yn (apply-contagion y precision))
  554 + (x (apply-contagion x precision))
  555 + xn z w a an pwr4 n epslon lambda sn s)
  556 + (cond ((and (zerop (imagpart yn))
  557 + (minusp (realpart yn)))
  558 + (setf xn (- x y))
  559 + (setf yn (- yn))
  560 + (setf z yn)
  561 + (setf w (sqrt (/ x xn))))
  562 + (t
  563 + (setf xn x)
  564 + (setf z yn)
  565 + (setf w 1)))
  566 + (setf a (/ (+ xn yn yn) 3))
  567 + (setf epslon (/ (abs (- a xn)) (errtol x y)))
  568 + (setf an a)
  569 + (setf pwr4 1)
  570 + (setf n 0)
  571 + (loop while (> (* epslon pwr4) (abs an))
  572 + do
  573 + (setf pwr4 (/ pwr4 4))
  574 + (setf lambda (+ (* 2 (sqrt xn) (sqrt yn)) yn))
  575 + (setf an (/ (+ an lambda) 4))
  576 + (setf xn (/ (+ xn lambda) 4))
  577 + (setf yn (/ (+ yn lambda) 4))
  578 + (incf n))
  579 + ;; c2=3/10,c3=1/7,c4=3/8,c5=9/22,c6=159/208,c7=9/8
  580 + (setf sn (/ (* pwr4 (- z a)) an))
  581 + (setf s (* sn sn (+ 3/10
  582 + (* sn (+ 1/7
  583 + (* sn (+ 3/8
  584 + (* sn (+ 9/22
  585 + (* sn (+ 159/208
  586 + (* sn 9/8))))))))))))
  587 + (/ (* w (+ 1 s))
  588 + (sqrt an))))
  589 +
  590 +(defun carlson-rj1 (x y z p)
  591 + (let* ((xn x)
  592 + (yn y)
  593 + (zn z)
  594 + (pn p)
  595 + (en (* (- pn xn)
  596 + (- pn yn)
  597 + (- pn zn)))
  598 + (sigma 0)
  599 + (power4 1)
  600 + (k 0)
  601 + (a (/ (+ xn yn zn pn pn) 5))
  602 + (epslon (/ (max (abs (- a xn))
  603 + (abs (- a yn))
  604 + (abs (- a zn))
  605 + (abs (- a pn)))
  606 + (errtol x y z p)))
  607 + (an a)
  608 + xnroot ynroot znroot pnroot lam dn)
  609 + (loop while (> (* power4 epslon) (abs an))
  610 + do
  611 + (setf xnroot (sqrt xn))
  612 + (setf ynroot (sqrt yn))
  613 + (setf znroot (sqrt zn))
  614 + (setf pnroot (sqrt pn))
  615 + (setf lam (+ (* xnroot ynroot)
  616 + (* xnroot znroot)
  617 + (* ynroot znroot)))
  618 + (setf dn (* (+ pnroot xnroot)
  619 + (+ pnroot ynroot)
  620 + (+ pnroot znroot)))
  621 + (setf sigma (+ sigma
  622 + (/ (* power4
  623 + (carlson-rc 1 (+ 1 (/ en (* dn dn)))))
  624 + dn)))
  625 + (setf power4 (* power4 1/4))
  626 + (setf en (/ en 64))
  627 + (setf xn (* (+ xn lam) 1/4))
  628 + (setf yn (* (+ yn lam) 1/4))
  629 + (setf zn (* (+ zn lam) 1/4))
  630 + (setf pn (* (+ pn lam) 1/4))
  631 + (setf an (* (+ an lam) 1/4))
  632 + (incf k))
  633 + (let* ((xndev (/ (* (- a x) power4) an))
  634 + (yndev (/ (* (- a y) power4) an))
  635 + (zndev (/ (* (- a z) power4) an))
  636 + (pndev (* -0.5 (+ xndev yndev zndev)))
  637 + (ee2 (+ (* xndev yndev)
  638 + (* xndev zndev)
  639 + (* yndev zndev)
  640 + (* -3 pndev pndev)))
  641 + (ee3 (+ (* xndev yndev zndev)
  642 + (* 2 ee2 pndev)
  643 + (* 4 pndev pndev pndev)))
  644 + (ee4 (* (+ (* 2 xndev yndev zndev)
  645 + (* ee2 pndev)
  646 + (* 3 pndev pndev pndev))
  647 + pndev))
  648 + (ee5 (* xndev yndev zndev pndev pndev))
  649 + (s (+ 1
  650 + (* -3/14 ee2)
  651 + (* 1/6 ee3)
  652 + (* 9/88 ee2 ee2)
  653 + (* -3/22 ee4)
  654 + (* -9/52 ee2 ee3)
  655 + (* 3/26 ee5)
  656 + (* -1/16 ee2 ee2 ee2)
  657 + (* 3/10 ee3 ee3)
  658 + (* 3/20 ee2 ee4)
  659 + (* 45/272 ee2 ee2 ee3)
  660 + (* -9/68 (+ (* ee2 ee5) (* ee3 ee4))))))
  661 + (+ (* 6 sigma)
  662 + (/ (* power4 s)
  663 + (sqrt (* an an an)))))))
  664 +
  665 +(defun carlson-rj (x y z p)
  666 + "Compute Carlson's Rj function:
  667 +
  668 + 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)"
  669 + (let* ((precision (float-contagion x y z p))
  670 + (xn (apply-contagion x precision))
  671 + (yn (apply-contagion y precision))
  672 + (zn (apply-contagion z precision))
  673 + (p (apply-contagion p precision))
  674 + (qn (- p)))
  675 + (cond ((and (and (zerop (imagpart xn)) (>= (realpart xn) 0))
  676 + (and (zerop (imagpart yn)) (>= (realpart yn) 0))
  677 + (and (zerop (imagpart zn)) (>= (realpart zn) 0))
  678 + (and (zerop (imagpart qn)) (> (realpart qn) 0)))
  679 + (destructuring-bind (xn yn zn)
  680 + (sort (list xn yn zn) #'<)
  681 + (let* ((pn (+ yn (* (- zn yn) (/ (- yn xn) (+ yn qn)))))
  682 + (s (- (* (- pn yn) (carlson-rj1 xn yn zn pn))
  683 + (* 3 (carlson-rf xn yn zn)))))
  684 + (setf s (+ s (* 3 (sqrt (/ (* xn yn zn)
  685 + (+ (* xn zn) (* pn qn))))
  686 + (carlson-rc (+ (* xn zn) (* pn qn)) (* pn qn)))))
  687 + (/ s (+ yn qn)))))
  688 + (t
  689 + (carlson-rj1 x y z p)))))
  690 +
  691 +;; Elliptic integral of the third kind:
  692 +;;
  693 +;; (A&S 17.2.14)
  694 +;;
  695 +;; PI(n; phi|m) = integrate(1/sqrt(1-m*sin(x)^2)/(1-n*sin(x)^2), x, 0, phi)
  696 +;;
  697 +(defun elliptic-pi (n phi m)
  698 + "Compute elliptic integral of the third kind:
  699 +
  700 + PI(n; phi|m) = integrate(1/sqrt(1-m*sin(x)^2)/(1-n*sin(x)^2), x, 0, phi)"
  701 + ;; Note: Carlson's DRJ has n defined as the negative of the n given
  702 + ;; in A&S.
  703 + (let* ((precision (float-contagion n phi m))
  704 + (n (apply-contagion n precision))
  705 + (phi (apply-contagion phi precision))
  706 + (m (apply-contagion m precision))
  707 + (nn (- n))
  708 + (sin-phi (sin phi))
  709 + (cos-phi (cos phi))
  710 + (k (sqrt m))
  711 + (k2sin (* (- 1 (* k sin-phi))
  712 + (+ 1 (* k sin-phi)))))
  713 + (- (* sin-phi (carlson-rf (expt cos-phi 2) k2sin 1))
  714 + (* (/ nn 3) (expt sin-phi 3)
  715 + (carlson-rj (expt cos-phi 2) k2sin 1
  716 + (+ 1 (* nn (expt sin-phi 2))))))))
... ...
rt-tests.lisp
... ... @@ -928,4 +928,143 @@
928 928 (let ((rf (carlson-rf 0 2 #q1q0))
929 929 (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0))
930 930 (check-accuracy 212 rf true))
931   - nil)
932 931 \ No newline at end of file
  932 + nil)
  933 +
  934 +;; Elliptic integral of the third kind
  935 +
  936 +;; elliptic-pi(0,phi,m) = elliptic-f(phi, m)
  937 +(rt:deftest oct.elliptic-pi.1d
  938 + (loop for k from 0 to 100
  939 + for phi = (random (/ pi 2))
  940 + for m = (random 1d0)
  941 + for epi = (elliptic-pi 0 phi m)
  942 + for ef = (elliptic-f phi m)
  943 + for result = (check-accuracy 53 epi ef)
  944 + unless (eq nil result)
  945 + append (list (list phi m) result))
  946 + nil)
  947 +
  948 +(rt:deftest oct.elliptic-pi.1q
  949 + (loop for k from 0 below 100
  950 + for phi = (random (/ +pi+ 2))
  951 + for m = (random #q1)
  952 + for epi = (elliptic-pi 0 phi m)
  953 + for ef = (elliptic-f phi m)
  954 + for result = (check-accuracy 53 epi ef)
  955 + unless (eq nil result)
  956 + append (list (list phi m) result))
  957 + nil)
  958 +
  959 +;; DLMF 19.6.3
  960 +;;
  961 +;; PI(n; pi/2 | 0) = pi/(2*sqrt(1-n))
  962 +(rt:deftest oct.elliptic-pi.19.6.3.d
  963 + (loop for k from 0 below 100
  964 + for n = (random 1d0)
  965 + for epi = (elliptic-pi n (/ pi 2) 0)
  966 + for true = (/ pi (* 2 (sqrt (- 1 n))))
  967 + for result = (check-accuracy 49 epi true)
  968 + unless (eq nil result)
  969 + append (list (list (list k n) result)))
  970 + nil)
  971 +
  972 +(rt:deftest oct.elliptic-pi.19.6.3.q
  973 + (loop for k from 0 below 100
  974 + for n = (random #q1)
  975 + for epi = (elliptic-pi n (/ (float-pi n) 2) 0)
  976 + for true = (/ (float-pi n) (* 2 (sqrt (- 1 n))))
  977 + for result = (check-accuracy 210 epi true)
  978 + unless (eq nil result)
  979 + append (list (list (list k n) result)))
  980 + nil)
  981 +
  982 +#+nil
  983 +(rt:deftest oct.elliptic-pi.19.6.2.d
  984 + (loop for k from 0 below 100
  985 + for n = (random 1d0)
  986 + for epi = (elliptic-pi (- n) (/ (float-pi n) 2) n)
  987 + for true = (+ (/ (float-pi n) 4 (sqrt (+ 1 (sqrt n))))
  988 + (/ (elliptic-k n) 2))
  989 + for result = (check-accuracy 53 epi true)
  990 + when result
  991 + append (list (list (list k n) result)))
  992 + nil)
  993 +
  994 +
  995 +#||
  996 +;; elliptic-pi(n, phi, 0) =
  997 +;; atanh(sqrt(1-n)*tan(phi))/sqrt(1-n) n < 1
  998 +;; atanh(sqrt(n-1)*tan(phi))/sqrt(n-1) n > 1
  999 +;; tan(phi) n = 1
  1000 +(rt:deftest oct.elliptic-pi.n0.d
  1001 + (loop for k from 0 below 100
  1002 + for phi = (random (/ pi 2))
  1003 + for n = (random 1d0)
  1004 + for epi = (elliptic-pi n phi 0)
  1005 + for true = (/ (atanh (* (tan phi) (sqrt (- 1 n))))
  1006 + (sqrt (- 1 n)))
  1007 + for result = (check-accuracy 53 epi true)
  1008 + unless (eq nil result)
  1009 + append (list (list (list k n phi) result)))
  1010 + nil)
  1011 +
  1012 +(rt:deftest oct.elliptic-pi.n1.d
  1013 + (loop for k from 0 below 100
  1014 + for phi = (random (/ pi 2))
  1015 + for epi = (elliptic-pi 0 phi 0)
  1016 + for true = (tan phi)
  1017 + for result = (check-accuracy 53 epi true)
  1018 + unless (eq nil result)
  1019 + append (list (list (list k phi) result)))
  1020 + nil)
  1021 +
  1022 +(rt:deftest oct.elliptic-pi.n2.d
  1023 + (loop for k from 0 below 100
  1024 + for phi = (random (/ pi 2))
  1025 + for n = (+ 1d0 (random 100d0))
  1026 + for epi = (elliptic-pi n phi 0)
  1027 + for true = (/ (atanh (* (tan phi) (sqrt (- n 1))))
  1028 + (sqrt (- n 1)))
  1029 + for result = (check-accuracy 52 epi true)
  1030 + ;; Not sure if this formula holds when atanh gives a complex
  1031 + ;; result. Wolfram doesn't say
  1032 + when (and (not (complexp true)) result)
  1033 + append (list (list (list k n phi) result)))
  1034 + nil)
  1035 +
  1036 +(rt:deftest oct.elliptic-pi.n0.q
  1037 + (loop for k from 0 below 100
  1038 + for phi = (random (/ +pi+ 2))
  1039 + for n = (random #q1)
  1040 + for epi = (elliptic-pi n phi 0)
  1041 + for true = (/ (atanh (* (tan phi) (sqrt (- 1 n))))
  1042 + (sqrt (- 1 n)))
  1043 + for result = (check-accuracy 212 epi true)
  1044 + unless (eq nil result)
  1045 + append (list (list (list k n phi) result)))
  1046 + nil)
  1047 +
  1048 +(rt:deftest oct.elliptic-pi.n1.q
  1049 + (loop for k from 0 below 100
  1050 + for phi = (random (/ +pi+ 2))
  1051 + for epi = (elliptic-pi 0 phi 0)
  1052 + for true = (tan phi)
  1053 + for result = (check-accuracy 212 epi true)
  1054 + unless (eq nil result)
  1055 + append (list (list (list k phi) result)))
  1056 + nil)
  1057 +
  1058 +(rt:deftest oct.elliptic-pi.n2.q
  1059 + (loop for k from 0 below 100
  1060 + for phi = (random (/ +pi+ 2))
  1061 + for n = (+ #q1 (random #q1))
  1062 + for epi = (elliptic-pi n phi 0)
  1063 + for true = (/ (atanh (* (tan phi) (sqrt (- n 1))))
  1064 + (sqrt (- n 1)))
  1065 + for result = (check-accuracy 209 epi true)
  1066 + ;; Not sure if this formula holds when atanh gives a complex
  1067 + ;; result. Wolfram doesn't say
  1068 + when (and (not (complexp true)) result)
  1069 + append (list (list (list k n phi) result)))
  1070 + nil)
  1071 +||#
... ...