Changeset 390a7483f5658fe802d5d239070cdaa573adf4a5
 Timestamp:
 03/12/11 10:35:18 (3 years ago)
 Children:
 5bd5df9360268fae90fc41fcdf86b728f8a54e86
 Parents:
 147fa2c7aa5e99988a7c3fb35800865d22efbc51
 gitcommitter:
 Raymond Toy <toy.raymond@…> (03/12/11 10:35:18)
 Files:

 2 modified
Legend:
 Unmodified
 Added
 Removed

qdelliptic.lisp
re2d8d6 r390a74 409 409 (/ (float +pi+ m) 2)) 410 410 (t 411 (let ((precision (floatcontagion m))) 412 (carlsonrf 0 ( 1 m) 1))))) 411 (carlsonrf 0 ( 1 m) 1)))) 413 412 414 413 ;; Elliptic integral of the first kind. This is computed using … … 533 532 (* (/ m 3) 534 533 (carlsonrd 0 y 1))))))) 534 535 ;; Carlson's Rc function. 536 ;; 537 ;; Some interesting identities: 538 ;; 539 ;; log(x) = (x1)*rc(((1+x)/2)^2, x), x > 0 540 ;; asin(x) = x * rc(1x^2, 1), x<= 1 541 ;; acos(x) = sqrt(1x^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^21) * rc(x^2,1), x >= 1 545 ;; atanh(x) = x * rc(1,1x^2), x<=1 546 ;; 547 548 (defun carlsonrc (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 (floatcontagion x y)) 553 (yn (applycontagion y precision)) 554 (x (applycontagion 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 carlsonrj1 (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 (carlsonrc 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 carlsonrj (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 (floatcontagion x y z p)) 670 (xn (applycontagion x precision)) 671 (yn (applycontagion y precision)) 672 (zn (applycontagion z precision)) 673 (p (applycontagion 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 (destructuringbind (xn yn zn) 680 (sort (list xn yn zn) #'<) 681 (let* ((pn (+ yn (* ( zn yn) (/ ( yn xn) (+ yn qn))))) 682 (s ( (* ( pn yn) (carlsonrj1 xn yn zn pn)) 683 (* 3 (carlsonrf xn yn zn))))) 684 (setf s (+ s (* 3 (sqrt (/ (* xn yn zn) 685 (+ (* xn zn) (* pn qn)))) 686 (carlsonrc (+ (* xn zn) (* pn qn)) (* pn qn))))) 687 (/ s (+ yn qn))))) 688 (t 689 (carlsonrj1 x y z p))))) 690 691 ;; Elliptic integral of the third kind: 692 ;; 693 ;; (A&S 17.2.14) 694 ;; 695 ;; PI(n; phim) = integrate(1/sqrt(1m*sin(x)^2)/(1n*sin(x)^2), x, 0, phi) 696 ;; 697 (defun ellipticpi (n phi m) 698 "Compute elliptic integral of the third kind: 699 700 PI(n; phim) = integrate(1/sqrt(1m*sin(x)^2)/(1n*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 (floatcontagion n phi m)) 704 (n (applycontagion n precision)) 705 (phi (applycontagion phi precision)) 706 (m (applycontagion m precision)) 707 (nn ( n)) 708 (sinphi (sin phi)) 709 (cosphi (cos phi)) 710 (k (sqrt m)) 711 (k2sin (* ( 1 (* k sinphi)) 712 (+ 1 (* k sinphi))))) 713 ( (* sinphi (carlsonrf (expt cosphi 2) k2sin 1)) 714 (* (/ nn 3) (expt sinphi 3) 715 (carlsonrj (expt cosphi 2) k2sin 1 716 (+ 1 (* nn (expt sinphi 2)))))))) 
rttests.lisp
r147fa2 r390a74 930 930 (checkaccuracy 212 rf true)) 931 931 nil) 932 933 ;; Elliptic integral of the third kind 934 935 ;; ellipticpi(0,phi,m) = ellipticf(phi, m) 936 (rt:deftest oct.ellipticpi.1d 937 (loop for k from 0 to 100 938 for phi = (random (/ pi 2)) 939 for m = (random 1d0) 940 for epi = (ellipticpi 0 phi m) 941 for ef = (ellipticf phi m) 942 for result = (checkaccuracy 53 epi ef) 943 unless (eq nil result) 944 append (list (list phi m) result)) 945 nil) 946 947 (rt:deftest oct.ellipticpi.1q 948 (loop for k from 0 below 100 949 for phi = (random (/ +pi+ 2)) 950 for m = (random #q1) 951 for epi = (ellipticpi 0 phi m) 952 for ef = (ellipticf phi m) 953 for result = (checkaccuracy 53 epi ef) 954 unless (eq nil result) 955 append (list (list phi m) result)) 956 nil) 957 958 ;; DLMF 19.6.3 959 ;; 960 ;; PI(n; pi/2  0) = pi/(2*sqrt(1n)) 961 (rt:deftest oct.ellipticpi.19.6.3.d 962 (loop for k from 0 below 100 963 for n = (random 1d0) 964 for epi = (ellipticpi n (/ pi 2) 0) 965 for true = (/ pi (* 2 (sqrt ( 1 n)))) 966 for result = (checkaccuracy 49 epi true) 967 unless (eq nil result) 968 append (list (list (list k n) result))) 969 nil) 970 971 (rt:deftest oct.ellipticpi.19.6.3.q 972 (loop for k from 0 below 100 973 for n = (random #q1) 974 for epi = (ellipticpi n (/ (floatpi n) 2) 0) 975 for true = (/ (floatpi n) (* 2 (sqrt ( 1 n)))) 976 for result = (checkaccuracy 210 epi true) 977 unless (eq nil result) 978 append (list (list (list k n) result))) 979 nil) 980 981 #+nil 982 (rt:deftest oct.ellipticpi.19.6.2.d 983 (loop for k from 0 below 100 984 for n = (random 1d0) 985 for epi = (ellipticpi ( n) (/ (floatpi n) 2) n) 986 for true = (+ (/ (floatpi n) 4 (sqrt (+ 1 (sqrt n)))) 987 (/ (elliptick n) 2)) 988 for result = (checkaccuracy 53 epi true) 989 when result 990 append (list (list (list k n) result))) 991 nil) 992 993 994 # 995 ;; ellipticpi(n, phi, 0) = 996 ;; atanh(sqrt(1n)*tan(phi))/sqrt(1n) n < 1 997 ;; atanh(sqrt(n1)*tan(phi))/sqrt(n1) n > 1 998 ;; tan(phi) n = 1 999 (rt:deftest oct.ellipticpi.n0.d 1000 (loop for k from 0 below 100 1001 for phi = (random (/ pi 2)) 1002 for n = (random 1d0) 1003 for epi = (ellipticpi n phi 0) 1004 for true = (/ (atanh (* (tan phi) (sqrt ( 1 n)))) 1005 (sqrt ( 1 n))) 1006 for result = (checkaccuracy 53 epi true) 1007 unless (eq nil result) 1008 append (list (list (list k n phi) result))) 1009 nil) 1010 1011 (rt:deftest oct.ellipticpi.n1.d 1012 (loop for k from 0 below 100 1013 for phi = (random (/ pi 2)) 1014 for epi = (ellipticpi 0 phi 0) 1015 for true = (tan phi) 1016 for result = (checkaccuracy 53 epi true) 1017 unless (eq nil result) 1018 append (list (list (list k phi) result))) 1019 nil) 1020 1021 (rt:deftest oct.ellipticpi.n2.d 1022 (loop for k from 0 below 100 1023 for phi = (random (/ pi 2)) 1024 for n = (+ 1d0 (random 100d0)) 1025 for epi = (ellipticpi n phi 0) 1026 for true = (/ (atanh (* (tan phi) (sqrt ( n 1)))) 1027 (sqrt ( n 1))) 1028 for result = (checkaccuracy 52 epi true) 1029 ;; Not sure if this formula holds when atanh gives a complex 1030 ;; result. Wolfram doesn't say 1031 when (and (not (complexp true)) result) 1032 append (list (list (list k n phi) result))) 1033 nil) 1034 1035 (rt:deftest oct.ellipticpi.n0.q 1036 (loop for k from 0 below 100 1037 for phi = (random (/ +pi+ 2)) 1038 for n = (random #q1) 1039 for epi = (ellipticpi n phi 0) 1040 for true = (/ (atanh (* (tan phi) (sqrt ( 1 n)))) 1041 (sqrt ( 1 n))) 1042 for result = (checkaccuracy 212 epi true) 1043 unless (eq nil result) 1044 append (list (list (list k n phi) result))) 1045 nil) 1046 1047 (rt:deftest oct.ellipticpi.n1.q 1048 (loop for k from 0 below 100 1049 for phi = (random (/ +pi+ 2)) 1050 for epi = (ellipticpi 0 phi 0) 1051 for true = (tan phi) 1052 for result = (checkaccuracy 212 epi true) 1053 unless (eq nil result) 1054 append (list (list (list k phi) result))) 1055 nil) 1056 1057 (rt:deftest oct.ellipticpi.n2.q 1058 (loop for k from 0 below 100 1059 for phi = (random (/ +pi+ 2)) 1060 for n = (+ #q1 (random #q1)) 1061 for epi = (ellipticpi n phi 0) 1062 for true = (/ (atanh (* (tan phi) (sqrt ( n 1)))) 1063 (sqrt ( n 1))) 1064 for result = (checkaccuracy 209 epi true) 1065 ;; Not sure if this formula holds when atanh gives a complex 1066 ;; result. Wolfram doesn't say 1067 when (and (not (complexp true)) result) 1068 append (list (list (list k n phi) result))) 1069 nil) 1070 #