modeq:=proc(lev,p)
local x2,x3,x4,x,alpha,beta,u,v,coe,suma,sb,
ecus,ecucoef,formu,sol,pol,dpol,pot,val,num,nu: global h,P:
x2:=q->64*q/(64*q+product((1+q^n)^(-24),n=1..300)):
x3:=q->27*q/(27*q+product((1+q^n+q^(2*n))^(-12),n=1..300)):
x4:=q->16*q*product(((1+q^(2*n))/(1+q^(2*n-1)))^8,n=1..300):
if lev=2 then x:=q->x2(q): val:=simplify((p+1)/4): h:=4/denom(val): dpol:=2*val: fi:
if lev=3 then x:=q->x3(q): val:=simplify((p+1)/3):
h:=6/denom(val): dpol:=numer(val): fi:
if lev=4 then x:=q->x4(q): val:=simplify((p+1)/8):
h:=8/denom(val): dpol:=numer(val): fi:
coe:=seq(seq(c[k-i,i],i=0..k),k=1..dpol):
num:=nops({coe}): print(level=lev, degree=p):
alpha:=q->x(q^(p)): beta:=q->x(q): sb:=dpol*(dpol+1):
u:=q->series((alpha(q)*beta(q))^(1/h),q,sb+30):
v:=q->series(((1-alpha(q))*(1-beta(q)))^(1/h),q,sb+30):
suma:=sum(sum(c[k-i,i]*u(q)^i*v(q)^(k-i),i=0..k),k=1..dpol):
ecus:=series(suma,q,sb)-1:
ecucoef:=coeffs(simplify(convert(ecus,polynom)),q):
formu:=sum(sum(c[k-i,i]*u^i*v^(k-i),i=0..k),k=1..dpol)-1:
sol:=solve({ecucoef},{coe}):
pol:=sort(subs(sol,formu),[u,v]):
P:=(uu,vv)->subs({u=uu,v=vv},pol):
print(u^h=alpha*beta,v^h=(1-alpha)*(1-beta)):
print(pol=0):
end: