load("c:\\dfw5\\math\\elliptbf2.mac")$

/* Use numerical checks from Numerical Algorithms 10(1995)13-26.
B.C. Carlson Numerical computation of real or complex elliptic integrals
page 22 3. Numerical checks */

mrf:bfloat([[1,2,0,1.3110287771461],
[%i,-%i,0,1.8540746773014],
[0.5,1,0,1.8540746773014],
[%i-1,%i,0,0.79612586584234-%i*(1.2138566698365)],
[2,3,4,0.58408284167715],[%i,-%i,2,1.0441445654064],
[%i-1,%i,1-%i,0.93912050218619-%i*(0.53296252018635)]]);

fpprec:20;
bferrtol:1.b-6;
maxn:16;

/* rf(x,y,z):=quad_qag 
(u/sqrt((u^2*(x+1)-2*u*x+x)*(u^2*(y+1)-2*u*y+y)*(u^2*(z+1)-2*u*z+z)),u,0,1,1);

rfi:rf(2,3,4);
rfi[1]-float(mrf[5][4]);

rf(x,y,z):=romberg 
(u/sqrt((u^2*(x+1)-2*u*x+x)*(u^2*(y+1)-2*u*y+y)*(u^2*(z+1)-2*u*z+z)),u,0,1);

rfi:rf(2,3,4);

rfi-float(mrf[5][4]); */

for n:1 thru 7 do (
rfi:expand(bfloat(bfrf(mrf[n][1],mrf[n][2],mrf[n][3])-mrf[n][4])),
display(n),display(rfi));

/* Check R_C. */

mrc:bfloat([[0,1/4,%pi],[9/4,2,log(2)],[0,%i,(1-%i)*(1.1107207345396)],
[-%i,%i,1.2260849569072-%i*(0.34471136988768)],[1/4,-2,log(2)/3],
[%i,-1,0.77778596920447+%i*(0.19832484993429)]]);

for n:1 thru 6 do (
rfi:expand(bfloat(bfrc(mrc[n][1],mrc[n][2])-mrc[n][3])),
display(n),display(rfi));

/* I can add some more exact solutions. */

mrc:bfloat([[0,1/4,%pi],[9/4,2,log(2)],[2,1,-log(sqrt(2)-1)],
[-%i,%i,-log(sqrt(2)-1)/2+
%pi/4-%i*(log(sqrt(2)-1)/2+%pi/4)],[1/4,-2,log(2)/3],
[%i,-1,sqrt(sqrt(2)/4-1/4)*atan(sqrt(sqrt(2)-1))-
sqrt(sqrt(2)/16+1/16)*log(-sqrt(2*sqrt(2)+2)+sqrt(2)+1)+
%i*(sqrt(sqrt(2)/4+1/4)*atan(sqrt(sqrt(2)-1))+sqrt(sqrt(2)/16-
1/16)*log(-sqrt(2*sqrt(2)+2)+sqrt(2)+1))],[0,1,%pi/2],
[%i,%i+1,%pi/4+%i*log(sqrt(2)-1)/2]]);

for n:1 thru 7 do (
rfi:expand(bfloat(bfrc(mrc[n][1],mrc[n][2])-mrc[n][3])),
display(n),display(rfi));

/* Check R_J. */

mrj:bfloat([[0,1,2,3,0.77688623778582],[2,3,4,5,0.14297579667157],
[2,3,4,-1+%i,0.13613945827771-%i*(0.38207561624427)],
[%i,-%i,0,2,1.6490011662711],[-1+%i,-1-%i,1,2,0.9414835884122],
[%i,-%i,0,1-%i,1.8260115229009+%i*(1.2290661908643)],
[-1+%i,-1-%i,1,-3+%i,-0.61127970812028-%i*(1.0684038390007)],
[-1+%i,-2-%i,-%i,-1+%i,1.8249027393704-%i*(1.2218475784827)],
[2,3,4,-0.5,0.24723819703052],[2,3,4,-5,-0.12711230042964]]);

for n:1 thru 10 do (
rfi:expand(bfloat(bfrj(mrj[n][1],mrj[n][2],mrj[n][3],mrj[n][4])-mrj[n][5]
)),
display(n),display(rfi));

/* page 23 Check R_D. */

mrd:bfloat([[0,2,1,1.7972103521034],[2,3,4,0.16510527294261],
[%i,-%i,2,0.6593385415422],[0,%i,-%i,1.270819627191+%i*(2.7811120159521)]
,
[0,%i-1,%i,-1.8577235439239-%i*(0.96193450888839)],
[-2-%i,-%i,-1+%i,1.8249027393704-%i*(1.2218475784827)]]);

for n:1 thru 6 do (
rfi:expand(bfloat(bfrd(mrd[n][1],mrd[n][2],mrd[n][3])-mrd[n][4])),
display(n),display(rfi));

/* Check R_G. */

mrg:bfloat([[0,16,16,%pi],[2,3,4,1.7255030280692],[0,%i,-%i,0.42360654239
699],
[0,%i-1,%i,0.44660591677018+%i*(0.70768352357515)],
[-%i,%i-1,%i,0.36023392184473+%i*(0.40348623401722)],
[0,0.0796,4,1.0284758090288]]);

for n:1 thru 6 do (
rfi:expand(bfloat(bfrg(mrg[n][1],mrg[n][2],mrg[n][3])-mrg[n][4])),
display(n),display(rfi));

/* Check the Legendre elliptic integrals using numerical integration. */

for m:1/8 thru 7/8 step 1/8 do(
	 for phi:%pi/16 thru %pi/2 step %pi/16 do(
      rfi:float(bfloat(bfelliptic_f(phi,m))-elliptic_f(phi,m)),
      display(m),display(phi),display(rfi)
));

for m:1/8 thru 7/8 step 1/8 do(
	 for phi:%pi/16 thru %pi/2 step %pi/16 do(
      rfi:float(bfloat(bfelliptic_e(phi,m))-elliptic_e(phi,m)),
      display(m),display(phi),display(rfi)
));

for n:1/8 thru 7/8 step 1/8 do(
	for m:1/8 thru 7/8 step 1/8 do(
		 for phi:%pi/16 thru %pi/2 step %pi/16 do(
    	  rfi:float(bfloat(bfelliptic_pi(n,phi,m))-elliptic_pi(n,phi,m)),
      	display(n),display(m),display(phi),display(rfi)
)));