"Carlson.mth Carlson elliptic integrals."

"These numerical algorithms are based on Numerical Algorithms 10(1995)13-26."

"B.C. Carlson Numerical computation of real or complex elliptic integrals"

"They have been updated to the algorithms in Journal of Computational"

"and Applied Mathematics 118 (2000) 71-85"

"Reduction theorems for elliptic integrands with the square root of"

"two quadratic factors"

"B.C. Carlson, James FitzSimons"

"Algorithm for R_F"

"page 16 2.3"

"page 16 2.4"

"Iterative definition of R_F"

;Sub(User')
R_FS9(e2,e3):=195*e2^8/32768-429*e2^7/59392+3003*e2^6*e3/63488+231*e2^6/25600~
-273*e2^5*e3^2/2048-77*e2^5*e3/1536-3*e2^5/256+3465*e2^4*e3^2/29696+315*e2^4*~
e3/5888+35*e2^4/2176-1155*e2^3*e3^3/7936-63*e2^3*e3^2/640-35*e2^3*e3/608-5*e2~
^3/208+105*e2^2*e3^4/1024+35*e2^2*e3^3/384+5*e2^2*e3^2/64+e2^2*e3/16+e2^2/24-~
315*e2*e3^4/7424-35*e2*e3^3/736-15*e2*e3^2/272-3*e2*e3/44-e2/10+63*e3^5/7936+~
7*e3^4/640+5*e3^3/304+3*e3^2/104+e3/14

R_F(xx,yy,zz,v,v1,q,a,n,03bb,x,y,z,e2,e3):=PROG(v:=[xx,yy,zz],q:=VECTOR(ABS(~
x),x,v),a:=MIN(q),q:=MAX(q),n:=CEILING((LN(IF(a>0,q/a,q))+LN(10)*PrecisionDig~
its/16)/LN(4))+4,LOOP(v1:=VECTOR(SQRT(i_),i_,v),03bb:=v1 SUB 1*v1 SUB 2+v1 S~
UB 2*v1 SUB 3+v1 SUB 3*v1 SUB 1,v:=VECTOR(i_+03bb,i_,v)/4,n:-1,IF(n=0,exit))~
,a:=AVERAGE(v),x:=1-v SUB 1/a,y:=1-v SUB 2/a,z:=-x-y,e2:=x*y-z^2,e3:=x*y*z,a:~
=1/SQRT(a),a*R_FS9(e2,e3)+a)

"Algorithm for R_C"

"page 17 2.10"

"Iterative definition of R_C"

"Add a special case when y is real and negative."

"page 18 2.14"

R_C_AUX2(x):=SUM((-x)^i_/(2*i_+1),i_,1,1+ABS(PrecisionDigits*LN(10)/(-LN(ABS(~
x)+0.01))))

R_C_AUX5(xx,yy,x,y,03bb):=PROG(x:=xx,y:=yy,LOOP(IF(ABS((y-x)/x)<1/64,RETURN(~
(R_C_AUX2((y-x)/x)+1)/SQRT(x))),03bb:=2*SQRT(x)*SQRT(y)+y,x:=(x+03bb)/4,y:=~
(y+03bb)/4))

R_C(x,y):=IF(IM(y)=0 AND RE(y)<0,IF(x=0,0,SQRT(x/(x-y))*R_C_AUX5(x-y,-y)),IF(~
x=0,pi/(2*SQRT(y)),R_C_AUX5(x,y)))

"Algorithm for R_J"

"page 18 2.17"

"page 19 2.21"

"Iterative definition of R_J"

S_M2(x):=IF(ABS(x)<1/64,SUM((-x/4)^i_*COMB(2*i_,i_)/(1-2*i_),i_,1,1+ABS(Preci~
sionDigits*LN(10)/(-LN(ABS(x)+0.01)))),SQRT(1+x)-1)

;Sub(User')
R_JS9(e2,e3,e4,e5):=3861*e2^8/229376-1287*e2^7/63488+273*e2^6*e3/2048-1287*e2~
^6*e4/10240+77*e2^6/3072-3861*e2^5*e3^2/10240-2079*e2^5*e3/14848+2079*e2^5*e4~
/15872-63*e2^5*e5/512-189*e2^5/5888+10395*e2^4*e3^2/31744-315*e2^4*e3*e4/512+~
297*e2^4*e3*e5/512+189*e2^4*e3/1280+297*e2^4*e4^2/1024-35*e2^4*e4/256+945*e2^~
4*e5/7424+105*e2^4/2432-105*e2^3*e3^3/256+297*e2^3*e3^2*e4/256-35*e2^3*e3^2/1~
28+945*e2^3*e3*e4/1856-945*e2^3*e3*e5/1984-5*e2^3*e3/32-945*e2^3*e4^2/3968+31~
5*e2^3*e4*e5/704+105*e2^3*e4/736-27*e2^3*e5^2/128-21*e2^3*e5/160-e2^3/16+297*~
e2^2*e3^4/1024+945*e2^2*e3^3/3712-2835*e2^2*e3^2*e4/3968+945*e2^2*e3^2*e5/140~
8+315*e2^2*e3^2/1472+945*e2^2*e3*e4^2/1408-81*e2^2*e3*e4*e5/64-63*e2^2*e3*e4/~
160+35*e2^2*e3*e5/96+45*e2^2*e3/272-27*e2^2*e4^3/128+35*e2^2*e4^2/192-315*e2^~
2*e4*e5/928-45*e2^2*e4/304+315*e2^2*e5^2/1984+15*e2^2*e5/112+9*e2^2/88-945*e2~
*e3^4/7936+315*e2*e3^3*e4/704-27*e2*e3^3*e5/64-21*e2*e3^3/160-81*e2*e3^2*e4^2~
/128+35*e2*e3^2*e4/96-315*e2*e3^2*e5/928-45*e2*e3^2/304-315*e2*e3*e4^2/928+31~
5*e2*e3*e4*e5/496+15*e2*e3*e4/56-105*e2*e3*e5^2/352-45*e2*e3*e5/184-9*e2*e3/5~
2+105*e2*e4^3/992-105*e2*e4^2*e5/352-45*e2*e4^2/368+9*e2*e4*e5^2/32+9*e2*e4*e~
5/40+3*e2*e4/20-5*e2*e5^2/48-9*e2*e5/68-3*e2/14+63*e3^5/2816-27*e3^4*e4/256+3~
5*e3^4/1152-105*e3^3*e4/928+105*e3^3*e5/992+5*e3^3/112+315*e3^2*e4^2/1984-105~
*e3^2*e4*e5/352-45*e3^2*e4/368+9*e3^2*e5^2/64+9*e3^2*e5/80+3*e3^2/40-35*e3*e4~
^3/352+9*e3*e4^2*e5/32+9*e3*e4^2/80-5*e3*e4*e5/24-9*e3*e4/68+45*e3*e5^2/464+9~
*e3*e5/76+e3/6+3*e4^4/128-5*e4^3/144+45*e4^2*e5/464+9*e4^2/152-45*e4*e5^2/496~
-3*e4*e5/28-3*e4/22+5*e5^3/176+9*e5^2/184+3*e5/26

R_J_AUX6(v,q,a,n,v1,x,y,z,p,03bb,03b1,sm,03b4,m,a1,e2,e3,e4,e5):=PROG(q:=V~
ECTOR(ABS(x),x,v),a:=MIN(q),q:=MAX(q),n:=CEILING((LN(IF(a>0,q/a,q))+LN(10)*Pr~
ecisionDigits/16)/LN(4))+4,x:=SQRT(v SUB 1),y:=SQRT(v SUB 2),z:=SQRT(v SUB 3)~
,p:=v SUB 4,03bb:=x*y+y*z+z*x,v1:=VECTOR(i_+03bb,i_,v)/4,03b1:=p*(x+y+z)+x~
*y*z,03b4:=PRODUCT(p-v SUB i_,i_,1,3),sm:=03b1+03b1*S_M2(03b4/03b1^2)/2,~
03b4:=03b4/4,m:=1/4,LOOP(x:=SQRT(v1 SUB 1),y:=SQRT(v1 SUB 2),z:=SQRT(v1 SUB~
 3),p:=SQRT(v1 SUB 4),03bb:=x*y+y*z+z*x,dm:=(p+x)*(p+y)*(p+z),rm:=2*sm+sm*S_~
M2(03b4/sm^2),sm:=(dm*rm-m*03b4)/(2*(dm+m*rm)),03b4:=03b4/4,m:=m/4,v1:=VE~
CTOR(i_+03bb,i_,v1)/4,n:-1,IF(n<=0,exit)),x:=v1 SUB 1,y:=v1 SUB 2,z:=v1 SUB ~
3,p:=v1 SUB 4,a:=(x+y+z+2*p)/5,x:=1-x/a,y:=1-y/a,z:=1-z/a,p:=(-x-y-z)/2,e2:=-~
3*p^2+z*(y+x)+x*y,e3:=2*p*e2+4*p^3+x*y*z,e4:=p*(p*e2+3*p^3+2*x*y*z),e5:=x*y*z~
*p^2,a:=1/a^(3/2),y:=03b4/sm^2,m*(a*R_JS9(e2,e3,e4,e5)+a)+3/sm*IF(ABS(y)<1/6~
4,R_C_AUX2(y)+1,R_C_AUX5(1,1+y)))

"Add a special case when p is real and negative."

"page 19 2.26"

SORT_RJ(v):=IF(RE(v SUB 1)>RE(v SUB 2),SORT_RJ([v SUB 2,v SUB 1,v SUB 3,v SUB~
 4]),IF(RE(v SUB 2)>RE(v SUB 3),SORT_RJ([v SUB 1,v SUB 3,v SUB 2,v SUB 4]),[v~
 SUB 1,v SUB 2,v SUB 3,-v SUB 4,(v SUB 1*(v SUB 2-v SUB 3)+v SUB 2*(v SUB 3-v~
 SUB 4))/(v SUB 2-v SUB 4)]))

R_J(x,y,z,p,v):=IF(IM(p)=0 AND RE(p)<0,PROG(v:=SORT_RJ([x,y,z,p]),((v SUB 5-v~
 SUB 2)*R_J_AUX6([v SUB 1,v SUB 2,v SUB 3,v SUB 5])-3*R_F(v SUB 1,v SUB 2,v S~
UB 3)+3*SQRT(v SUB 1*v SUB 2*v SUB 3/(v SUB 1*v SUB 3+v SUB 5*v SUB 4))*R_C(v~
 SUB 1*v SUB 3+v SUB 5*v SUB 4,v SUB 5*v SUB 4))/(v SUB 2+v SUB 4)),R_J_AUX6(~
[x,y,z,p]))

"Algorithm for R_D"

"page 20 2.28"

"page 20 2.27"

"Iterative definition of R_D"

R_D(xx,yy,zz,v,v1,v2,q,a,x,y,z,03bb,s,f,e2,e3,e4,e5):=PROG(v:=[xx,yy,zz],q:=~
VECTOR(ABS(x),x,v),a:=MIN(q),q:=MAX(q),n:=CEILING((LN(IF(a>0,q/a,q))+LN(10)*P~
recisionDigits/16)/LN(4))+4,s:=0,f:=1,v1:=v,LOOP(x:=SQRT(v1 SUB 1),y:=SQRT(v1~
 SUB 2),z:=SQRT(v1 SUB 3),03bb:=x*y+y*z+z*x,s:+f*3/(z*(v1 SUB 3+03bb)),v1:=~
VECTOR(i_+03bb,i_,v1)/4,f:/4,n:-1,IF(n=0,exit)),x:=v1 SUB 1,y:=v1 SUB 2,z:=v~
1 SUB 3,a:=(x+y+3*z)/5,x:=1-x/a,y:=1-y/a,z:=(-x-y)/3,e2:=x*y-6*z^2,e3:=(3*x*y~
-8*z^2)*z,e4:=3*(x*y-z^2)*z^2,e5:=x*y*z^3,a:=1/a^(3/2),f*(a*R_JS9(e2,e3,e4,e5~
)+a)+s)

"Algorithm for R_G"

"page 15 1.7"

R_G(x,y,z):=(z*R_F(x,y,z)-(x-z)*(y-z)*R_D(x,y,z)/3+SQRT(x*y/z))/2

"page 24 4. Other integrals"

"These are the Legendre elliptic integrals."

LEGENDRE_F(03c6,m):=SIN(03c6)*R_F(COS(03c6)^2,1-m*SIN(03c6)^2,1)

LEGENDRE_E1(03c6,m):=LEGENDRE_F(03c6,m)-m/3*SIN(03c6)^3*R_D(COS(03c6)^2,1~
-m*SIN(03c6)^2,1)

LEGENDRE_E(03c6,m):=IF(0<=m<=1 AND m*SIN(03c6)^2/=1,(1-m)*LEGENDRE_F(03c6,~
m)+m*(1-m)/3*SIN(03c6)^3*R_D(COS(03c6)^2,1,1-m*SIN(03c6)^2)+m*SIN(03c6)*C~
OS(03c6)/SQRT(1-m*SIN(03c6)^2),IF(1<=m<=1/SIN(03c6)^2 AND 03c6/=pi/2,(m-1~
)/3*SIN(03c6)^3*R_D(1-m*SIN(03c6)^2,1,COS(03c6)^2)+TAN(03c6)*SQRT(1-m*SIN~
(03c6)^2),LEGENDRE_F(03c6,m)-m/3*SIN(03c6)^3*R_D(COS(03c6)^2,1-m*SIN(03c~
6)^2,1)))

LEGENDRE_PI(03c6,m,n):=LEGENDRE_F(03c6,m)+n/3*SIN(03c6)^3*R_J(COS(03c6)^2~
,1-m*SIN(03c6)^2,1,1-n*SIN(03c6)^2)

ELLIPTICF(x,k):=x*R_F(1-x^2,1-k^2*x^2,1)

ELLIPTICE(x,k):=ELLIPTICF(x,k)-k^2*x^3/3*R_D(1-x^2,1-k^2*x^2,1)

ELLIPTICPI(x,n,k):=ELLIPTICF(x,k)+n*x^3/3*R_J(1-x^2,1-k^2*x^2,1,1-n*x^2)

ELLIPTICK(k):=R_F(0,1-k^2,1)

ELLIPTICEK(k):=ELLIPTICK(k)-k^2/3*R_D(0,1-k^2,1)

ELLIPTICPIK(n,k):=ELLIPTICK(k)+n/3*R_J(0,1-k^2,1,1-n)