/* elliptbf.mac jim fitzsimons 25 january, 1998
   for complex numbers */

bfrc(x, y):= block([xn,yn,z,w,a,an,pwr4,n,epslon,lamda,sn,s],
yn:expand(bfloat(y)),
if ((imagpart(yn) = 0) and (realpart(yn) < 0.b0)) then
   (xn:bfloat(x-y),
   yn:-yn,
   z:yn,
   w:sqrt(bfloat(x)/xn))
else
   (xn:bfloat(x),
   z:yn,
   w:1.b0),

a:bfloat((xn+yn+yn)/3.b0),
epslon:cabs(a-xn)/bferrtol,
an:a,
pwr4:1.b0,
n:0,
while (epslon*pwr4 > cabs(an)) do block([],
   pwr4:pwr4/4.b0,
   lamda:expand(bfloat(2.b0*sqrt(xn)*sqrt(yn)+yn)),
   an:expand(bfloat((an+lamda)/4.b0)),
   xn:expand(bfloat((xn+lamda)/4.b0)),
   yn:expand(bfloat((yn+lamda)/4.b0)),
   n:n+1),
/* c2=3/10,c3=1/7,c4=3/8,c5=9/22,c6=159/208,c7=9/8 */
sn:(z-a)*pwr4/an,
s:sn*sn*(3.b0/10.b0+sn*(1.b0/7.b0+sn*(3.b0/8.b0+
sn*(9.b0/22.b0+sn*(159.b0/208.b0+9.b0*sn/8.b0))))),
rectform(expand(bfloat(w*(1.b0+s)/sqrt(an)))) )$

bfrd(x, y, z):= block([xn,yn,zn,a,epslon,an,sigma,power4,n,
 xnroot,ynroot,znroot,lamda,xndev,yndev,zndev,ee2,ee3,ee4,ee5,s],
xn:bfloat(x),
yn:bfloat(y),
zn:bfloat(z),
a:bfloat((xn+yn+3.b0*zn)/5.b0),
epslon:max(cabs(a-xn),cabs(a-yn),cabs(a-zn))/bferrtol,
an:a,
sigma:0.b0,
power4:1.b0,
n:0,
while (power4*epslon > cabs(an)) do block([],
   xnroot:bfloat(sqrt(xn)),
   ynroot:bfloat(sqrt(yn)),
   znroot:bfloat(sqrt(zn)),
   lamda:expand(xnroot*ynroot+xnroot*znroot+ynroot*znroot),
   sigma:expand(bfloat(sigma+power4/(znroot*(zn+lamda)))),
   power4:power4/4.b0,
   xn:expand(bfloat((xn+lamda)/4.b0)),
   yn:expand(bfloat((yn+lamda)/4.b0)),
   zn:expand(bfloat((zn+lamda)/4.b0)),
   an:expand(bfloat((an+lamda)/4.b0)),
   n:n+1),
/* 
c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26 
*/
xndev:(a-bfloat(x))*power4/an,
yndev:(a-bfloat(y))*power4/an,
zndev:(-xndev-yndev)/3.b0,
ee2:xndev*yndev-6.b0*zndev*zndev,
ee3:(3.b0*xndev*yndev-8.b0*zndev*zndev)*zndev,
ee4:3.b0*(xndev*yndev-zndev*zndev)*zndev*zndev,
ee5:xndev*yndev*zndev*zndev*zndev,
/* s:ee2*(-3.b0/14.b0+9.b0*ee2/88.b0-9.b0*ee3/52.b0),
s2:ee3/6.b0-3.b0*ee4/22.b0+3.b0*ee5/26.b0, */
s:1.b0-3.b0*ee2/14.b0+ee3/6.b0+9.b0*ee2*ee2/88.b0-3.b0*ee4/22.b0
-9.b0*ee2*ee3/52.b0+3.b0*ee5/26.b0
-1.b0*ee2*ee2*ee2/16.b0+3.b0*ee3*ee3/10.b0+3.b0*ee2*ee4/20.b0
+45.b0*ee2*ee2*ee3/272.b0-9.b0*(ee2*ee5+ee3*ee4)/68.b0,
rectform(expand(bfloat(3.b0*sigma+power4*s/(an^(3/2))))) )$

bfrf(x, y, z):= block([xn,yn,zn,a,epslon,an,pwr4,xnroot,ynroot,znroot,
lamda,xndev,yndev,zndev,ee2,ee3,s,n],
xn:bfloat(x),
yn:bfloat(y),
zn:bfloat(z),
a:bfloat((xn+yn+zn)/3.b0),
epslon:max(cabs(a-xn),cabs(a-yn),cabs(a-zn))/bferrtol,
an:a,
pwr4:1.b0,
n:0,
while ((epslon*pwr4 > cabs(an)) and (n < maxn)) do block([],
   xnroot:bfloat(sqrt(xn)),
   ynroot:bfloat(sqrt(yn)),
   znroot:bfloat(sqrt(zn)),
   lamda:expand(xnroot*ynroot+xnroot*znroot+ynroot*znroot),
   xn:expand(bfloat((xn+lamda)/4.b0)),
   yn:expand(bfloat((yn+lamda)/4.b0)),
   zn:expand(bfloat((zn+lamda)/4.b0)),
   an:expand(bfloat((an+lamda)/4.b0)),
   pwr4:pwr4/4.b0,
   n:n+1),
/* c1=-1/10,c2=1/14,c3=1/24,c4=-3/44 */
xndev:(a-bfloat(x))*pwr4/an,
yndev:(a-bfloat(y))*pwr4/an,
zndev:-xndev-yndev,
ee2:xndev*yndev-zndev*zndev,
ee3:xndev*yndev*zndev,
s:1.b0-1.b0*ee2/10.b0+ee3/14.b0+ee2*ee2/24.b0-3.b0*ee2*ee3/44.b0,
rectform(expand(bfloat(s/sqrt(an)))) )$

bfrj(x, y, z, p):= block([xn,yn,zn,pn,v,s],
xn:bfloat(x),
yn:bfloat(y),
zn:bfloat(z),
qn:bfloat(-p),
if ((imagpart(xn) = 0) and (realpart(xn) >= 0.b0) and
    (imagpart(yn) = 0) and (realpart(yn) >= 0.b0) and
    (imagpart(zn) = 0) and (realpart(zn) >= 0.b0) and
    (imagpart(qn) = 0) and (realpart(qn) > 0.b0)) then
   (v:sort([xn,yn,zn]),
    xn:v[1],yn:v[2],zn:v[3],
    pn:yn+(zn-yn)*(yn-xn)/(yn+qn),
    s:(pn-yn)*bfrj1(xn,yn,zn,pn)-3.b0*bfrf(xn,yn,zn),
    s:s+3.b0*sqrt(xn*yn*zn/(xn*zn+pn*qn))*bfrc(xn*zn+pn*qn,pn*qn),   
    rectform(expand(bfloat(s/(yn+qn)))) )   
else
   bfrj1(x, y, z, p))$

bfrj1(x, y, z, p):= block([xn,yn,zn,pn,sigma,power4,k,a,epslon,an,
xnroot,ynroot,znroot,pnroot,lamda,dn,en,xndev,yndev,zndev,pndev,
ee2,ee3,ee4,ee5,s],
xn:bfloat(x),
yn:bfloat(y),
zn:bfloat(z),
pn:bfloat(p),
en:expand(bfloat((pn-xn)*(pn-yn)*(pn-zn))),
sigma:0.b0,
power4:1.b0,
k:0,
a:bfloat((xn+yn+zn+pn+pn)/5.b0),
epslon:max(cabs(a-xn),cabs(a-yn),cabs(a-zn),cabs(a-pn))/bferrtol,
an:a,
while (power4*epslon > cabs(an)) do block([],
   xnroot:bfloat(sqrt(xn)),
   ynroot:bfloat(sqrt(yn)),
   znroot:bfloat(sqrt(zn)),
   pnroot:bfloat(sqrt(pn)),
   lamda:expand(xnroot*ynroot+xnroot*znroot+ynroot*znroot),
   dn:expand(bfloat((pnroot+xnroot)*(pnroot+ynroot)*(pnroot+znroot))),
   sigma:expand(bfloat(sigma+power4*bfrc(1.b0,1.b0+en/(dn*dn))/dn)),
   power4:power4/4.b0,
   en:expand(en/64.b0),
   xn:expand(bfloat((xn+lamda)/4.b0)),
   yn:expand(bfloat((yn+lamda)/4.b0)),
   zn:expand(bfloat((zn+lamda)/4.b0)),
   pn:expand(bfloat((pn+lamda)/4.b0)),
   an:expand(bfloat((an+lamda)/4.b0)),
   k:k+1),
/*1-3/14*ee2+1/6*ee3+9/88*ee2^2-3/22*ee4-9/52*ee2*ee3+3/26*e5
-1/16*ee2^3+3/10*ee3^2+3/20*ee2*ee4+45/272*ee2^2*ee3-9/68*(ee2*ee5+ee3*ee
4) */
xndev:(a-bfloat(x))*power4/an,
yndev:(a-bfloat(y))*power4/an,
zndev:(a-bfloat(z))*power4/an,
pndev:(-xndev-yndev-zndev)/2.b0,
ee2:xndev*yndev+xndev*zndev+yndev*zndev-3.b0*pndev*pndev,
ee3:xndev*yndev*zndev+2.b0*ee2*pndev+4.b0*pndev*pndev*pndev,
ee4:(2.b0*xndev*yndev*zndev+ee2*pndev+3.b0*pndev*pndev*pndev)*pndev,
ee5:xndev*yndev*zndev*pndev*pndev,
/* s:1.b0-3.b0*ee2/14.b0+ee3/6.b0+9.b0*ee2*ee2/88.b0-9.b0*ee2*ee3/52.b0
+3.b0*ee5/26.b0, */ 
s:1.b0-3.b0*ee2/14.b0+ee3/6.b0+9.b0*ee2*ee2/88.b0-3.b0*ee4/22.b0
-9.b0*ee2*ee3/52.b0+3.b0*ee5/26.b0
-1.b0*ee2*ee2*ee2/16.b0+3.b0*ee3*ee3/10.b0+3.b0*ee2*ee4/20.b0
+45.b0*ee2*ee2*ee3/272.b0-9.b0*(ee2*ee5+ee3*ee4)/68.b0,
rectform(expand(bfloat(6.b0*sigma+power4*s/(an^(3/2))))) )$

bfellf(phi, ak):=block([s,c,cc,sak,q],
s:sin(phi),
c:cos(phi),
cc:c * c,
sak:s * bfloat(ak),
q:(1.b0 - sak) * (1.b0 + sak),
rectform(bfloat(expand(s * bfrf(cc, q, 1.b0)))) )$

bfelle(phi, ak):=block([s,c,cc,sak,q],
s:sin(phi),
c:cos(phi),
cc:c * c,
sak:s * bfloat(ak),
q:(1.b0 - sak) * (1.b0 + sak),
rectform(expand(bfloat(s*(bfrf(cc, q, 1.b0)-sak*sak*bfrd(cc, q, 
1.b0)/3.b0)))) )$

bfellpi(phi, en, ak):=block([s,c,cc,sak,q,enss],
s:sin(phi),
enss:en * s * s,
c:cos(phi),
cc:c * c,
sak:s * bfloat(ak),
q:(1.b0 - sak) * (1.b0 + sak),
rectform(expand(bfloat(s*(bfrf(cc, q, 1.b0)-enss*bfrj(cc, q, 1.b0, 1.b0 
+ enss) 
/ 3.b0)))) )$

bfam(u, m):=block([a,b,c,k:bfloat(sqrt(m)),fac,phi,i,j,dela,
ee2:bferrtol*bferrtol],
array([a,b,c],256),
a[0]:1.b0,
b[0]:bfloat(sqrt((1.b0-k)*(1.b0+k))),
c[0]:k,
if (1.b0-k) > ee2 then
   (fac:1.b0,dela:1.b0,i:1,
   while dela > ee2 do block([],
      fac:fac*2.b0,
      a[i]:(a[i-1]+b[i-1])/2.b0,
      b[i]:bfloat(sqrt(a[i-1]*b[i-1])),
      c[i]:(a[i-1]-b[i-1])/2.b0,
      dela:abs(a[i]-a[i-1]),
      i:i+1),
   phi:u*fac*a[i-1],
   for j:i-1 step -1 thru 1 do
      (phi:(phi+asin(c[j]/a[j]*sin(phi)))/2.b0))
else
   (phi:2.b0*atan(exp(u))-bfloat(%pi/2.b0)),
rectform(expand(bfloat(phi))) )$

bferrtol:1.b-6;

maxn:16;