fullsimplify:=proc(expression)
combine(evalc(simplify(expand(combine(rationalize(radnormal(
expand(simplify(rationalize(simplify(
combine(radnormal(expand(expression)),radicals))))))),radicals)))),radicals):
end:

dat5:=proc()
global Q,R,P,H,u0,v0,dec,rulepol: dec:=5:
Q:=(u,v)->u^3+v^3:
R:=(u,v)->-u^2*v^2+4:
P:=(u,v)->Q(u,v)^2-u*v*R(u,v)^2:
H:=(1+3*sqrt(57))^(1/3):
u0:=(-1/192+sqrt(57)/64)*H^2-1/3*H+2/3:
v0:=((-sqrt(57)/128+1/384)*H^2+1/6*H+2/3)-
((sqrt(3)/384-sqrt(171)/128)*H^2-sqrt(3)/6*H)*I:
rulepol:={u^3=2*u^2-4*u+2,v(u)^3=2*v(u)^2-4*v(u)+2}:
end:

dat11:=proc()
global Q,R,P,H,u0,v0,dec,rulepol: dec:=11:
Q:=(u,v)->u^6+v^6:
R:=(u,v)->u^5*v^5-11*u^4*v^4+44*u^3*v^3-88*u^2*v^2+88*u*v-32:
P:=(u,v)->Q(u,v)^2-u*v*R(u,v)^2:
H:=(35+3*sqrt(129))^(1/3):
u0:=(-35/48+sqrt(129)/16)*H^2-1/3*H+4/3:
v0:=((35/96-sqrt(129)/32)*H^2+1/6*H+4/3)+
((-35*sqrt(3)/96+3*sqrt(43)/32)*H^2+sqrt(3)/6*H)*I:
rulepol:={u^3=4*u^2-4*u-2,v(u)^3=4*v(u)^2-4*v(u)-2}:
end:

dat17:=proc() global Q,R,P,H,u0,v0,dec,rulepol: dec:=17:
Q:=(u,v)->(u^9+v^9)+17*(u^8*v^5+u^5*v^8)+
119*(u^6*v^3+u^3*v^6)+272*(u^4*v+u*v^4):
R:=(u,v)->-(u^8*v^8)-34*(u^7*v+u*v^7)+
34*(u^6*v^6)+340*(u^4*v^4)+544*(u^2*v^2)-256:
P:=(u,v)->Q(u,v)^2-(u*v)*R(u,v)^2:
H:=(91+9*sqrt(201))^(1/3):
u0:=(91/1200-3*sqrt(201)/400)*H^2+1/3*H-2/3:
v0:=((3*sqrt(201)/800-91/2400)*H^2-1/6*H-2/3)+
((-9*sqrt(67)/800+91*sqrt(3)/2400)*H^2-sqrt(3)/6*H)*I:
rulepol:={u^3=-2*u^2-8*u+2,v(u)^3=-2*v(u)^2-8*v(u)+2}:
end:

dat41:=proc() global Q,R,P,H,u0,v0,dec,rulepol: dec:=41:
Q:=(u,v)->u^(21)+v^(21)+943*(u^(20)*v^(5)+u^(5)*v^(20))+
123*(u^(20)*v^(17)+u^(17)*v^(20))+
40713*(u^(19)*v^(10)+u^(10)*v^(19))+72939*(u^18*v^3+u^3*v^18)+
3772*(u^18*v^15+u^15*v^18)+
733531*(u^17*v^8+u^8*v^17)+15088*(u^16*v+u*v^16)+
339111*(u^16*v^13+u^13*v^16)-
6494359*(u^15*v^6+u^6*v^15)+3112310*(u^14*v^11+u^11*v^14)+
11736496*(u^13*v^4+u^4*v^13)
-36004437*(u^12*v^9+u^9*v^12)+10422528*(u^11*v^2+u^2*v^11)+
49796960*(u^10*v^7+u^7*v^10)
+86812416*(u^8*v^5+u^5*v^8)+15450112*(u^6*v^3+u^3*v^6)+
8060928*(u^4*v+u*v^4):
R:=(u,v)->41*(u^(20)*v^(8)+u^(8)*v^(20))-
u^(20)*v^(20)+574*(u^(19)*v+u*v^(19))-
4059*(u^(19)*v^(13)+u^(13)*v^(19))-155554*(u^(18)*v^(6)+
u^(6)*v^(18))+574*u^(18)*v^(18)-
160310*(u^17*v^11+u^11*v^17)+701100*(u^16*v^4+u^4*v^16)-
2050*u^16*v^16-1753160*(u^15*v^9+u^9*v^15)-
2488864*(u^14*v^2+u^2*v^14)-462726*u^14*v^14+
20156994*(u^13*v^7+u^7*v^13)+10496*(u^12+v^12)-
3571756*u^12*v^12-28050560*(u^11*v^5+u^5*v^11)+
45567400*u^10*v^10-41039360*(u^9*v^3+u^3*v^9)-
57148096*u^8*v^8-16625664*(u^7*v+u*v^7)-118457856*u^6*v^6-
8396800*u^4*v^4+37617664*u^2*v^2-1048576:
P:=(u,v)->Q(u,v)^2-(u*v)*R(u,v)^2:
H:=(467+33*sqrt(489))^(1/3):
u0:=(-467/13872+11*sqrt(489)/4624)*H^2-1/3*H+4/3:
v0:=((467/27744-11*sqrt(489)/9248)*H^2+1/6*H+4/3)+
((467*sqrt(3)/27744-33*sqrt(163)/9248)*H^2
-sqrt(3)/6*H)*I:
rulepol:={u^3=4*u^2-28*u-2,v(u)^3=4*v(u)^2-28*v(u)-2}:
end:

chud:=proc(dd)
local T,w0,dv0,alpha0,beta0,z0,dalpha0,dbeta0,m0,ddv0,ddalpha0,
ddbeta0,b0,a0,J0, atu0,w,y,y0,dv,ddv,dw,ddw,dy,ddy,a1,a2,o: o:=time():
if dd=5 then dat5() fi: if dd=11 then dat11() fi:
if dd=17 then dat17() fi: if dd=41 then dat41() fi:
atu0:={u=u0,v(u)=v0,diff(v(u),u)=dv0,diff(diff(v(u),u),u)=ddv0}:
print(Modified-Weber-Polynomial=P(u,v)):
T:=u->P(u,v(u)):
print(u[0]=u0); print(v[0]=v0):
dv:=solve(diff(T(u),u)=0,diff(v(u),u)): ddv:=diff(dv,u):
w:=u^(12)/(u^(12)-16)^3: y:=v(u)^(12)/(v(u)^(12)-16)^3:
dw:=diff(w,u): ddw:=diff(dw,u):
dy:=subs(diff(v(u),u)=dv,diff(y,u)):
ddy:=subs(diff(v(u),u)=dv,diff(dy,u)):
w0:=fullsimplify(subs(atu0,simplify(w,rulepol)));
simplify(expand(P(u0,v0))):
T:=u->P(u,v(u)): print(): print():
alpha0:=solve(al*(1-al)=simplify(432*w0),al)[1]:
beta0:=1-alpha0; print(alpha[0]=alpha0): print(beta[0]=beta0):
z0:=expand(4*alpha0*beta0): print(z[0]=z0): J0:=12^3/z0:
b0:=fullsimplify(sqrt(4*dec-1)*(1-2*alpha0)):
print(b[0]=b0): print(): print():
dv0:=fullsimplify(subs(atu0,simplify(dv,rulepol))):
print(diff(v(u),u)(u[0])=dv0):
dalpha0:=fullsimplify(432*subs(atu0,simplify(dw,rulepol))/(1-2*alpha0)):
print(diff(alpha(u),u)(u[0])=dalpha0):
dbeta0:=fullsimplify(432*subs(atu0,simplify(dy,rulepol))/(1-2*beta0)):
print(diff(beta(u),u)(u[0])=dbeta0):
m0:=fullsimplify(sqrt(fullsimplify(1/dec*dalpha0/dbeta0))):
print(m[0]=m0):
ddv0:=fullsimplify(subs(atu0,simplify(ddv,rulepol))):
print(diff(diff(v(u),u),u)(u[0])=ddv0):
ddalpha0:=fullsimplify(2*dalpha0^2+
432*subs(atu0,simplify(ddw,rulepol)))/(1-2*alpha0):
print(diff(diff(alpha(u),u),u)(u[0])=ddalpha0):
ddbeta0:=fullsimplify(2*dbeta0^2+432*fullsimplify(
subs(atu0,simplify(ddy,rulepol))))/(1-2*beta0):
print(diff(diff(beta(u),u),u)(u[0])=ddbeta0):
a1:=fullsimplify(m0/2*(ddalpha0/dalpha0^2-ddbeta0/(dalpha0*dbeta0))):
a2:=fullsimplify(1/2*(m0+1/(dec*m0))*(alpha0-beta0)/(alpha0*beta0)):
a0:=expand(-2*alpha0*beta0*dec*(a1+a2)):
print(z[0]=z0,b[0]=b0,a[0]=a0,J[0]=J0):
print(Sum((6*n)!/((3*n)!*(n!)^3)*(b0*n+a0)/J0^n,n=0..infinity)=1/Pi):
print(seconds=time()-o):
end: