Changeset 390a74
 Timestamp:
 03/12/11 18:35:18 (4 years ago)
 Branches:
 master
 Children:
 5bd5df
 Parents:
 147fa2
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

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