ACC := 140;
t := series(q^(-5) * mul( ( (1-q^n)/(1-q^(11*n)) )^12, n=1..ACC), q=0, ACC);
E4 := 1 + 240*add( n^3*q^n/(1-q^n), n=1..ACC):
Ds := q * mul( (1-q^n)^(24), n=1..ACC):
f := series(t^3 * E4^3 / Ds, q=0, ACC);
h := series(q*t*mul(1-q^(11*k), k=1..40) * add(combinat[numbpart](11*k+6)*q^k, k=0..40), q=0, 40);


# Compute the algebraic relation p in Q[x,y] between t and f using an Ansatz:
n := 5;
m := 16;
p := add(add( c[i,j]*x^i*y^j, i=0..min(m, max(m,n)-j)), j=0..n);
IsZero := eval(p, {x=t, y=f}):
SolveCoeffsZero := proc(IsZero, q, ACC)
	solve({coeffs(convert(series(IsZero, q=0, ACC),polynom),q)})
end:
p := eval(p, SolveCoeffsZero(IsZero, q, ACC)):
p := sort(collect(primpart(p,y),y, factor), y);


# Compute a basis v of all integral functions i for which i/x^d is
# integral at infinity. Then write h as a linear combination of v:
d := 1;
read NormalBasis:
v := UpToDegree(p, x, y, d):
H := add(c[i]*v[i], i=1..nops(v)):
IsZero := h - eval(H, {x=t,y=f}):
h := factor(eval(H, SolveCoeffsZero(IsZero, q, 40)));


# Shorten h with "Rational Univariate Representation"
ShortenAlg := proc(h, f,x,y) local d, RUR;
	d := diff(f,y);
	RUR := rem(h*d, f, y)/d;
	map(factor, convert(RUR, parfrac, y))
end:
h := ShortenAlg(h, p,x,y);


with(algcurves):
genus(p,x,y);


# The genus is 1. So the integral elements can have poles at
# infinity of order 0, 2, 3, 4, .... (the only gap is 1).
# Call those functions 1, z[2], z[3], ... then we should
# be able to write h as a polynomial in z[2], z[3].

# Compute the q-expansions of a basis v of some subspace of Q(x)[y]/(f).
# Adjust this basis to bring these q-expansions into Reduced Echelon Form.
# Return the adjusted basis, as well as q-expansions up to O(q^1).
qEchelonForm := proc(v,q, x,y, t,f, ACC)
	local n,m,w,i,j,M;
	n := nops(v);
	w := [seq(convert(series(eval(i, {x=t, y=f}), q=0, ACC),polynom), i=v)];
	m := min(map(ldegree,w,q));
	M := Matrix(n,1-m,[seq(seq(coeff(i,q,j),j=m..0),i=w)]);
	M := < M | LinearAlgebra:-IdentityMatrix(n) > ;
	M := LinearAlgebra:-ReducedRowEchelonForm( M );
	[seq([factor(add(M[i,1-m+j]*v[j], j=1..n)),
	 add(M[i,j-m+1]*q^j,j=m..0)], i=1..n)]
end:

v := qEchelonForm(v, q, x,y, t,f, 40):
for i in v do
	PoleOrder := -ldegree(i[2],q);
	z[PoleOrder] := ShortenAlg(i[1], p,x,y);
	zq[PoleOrder] := i[2] + O(q);
od;

# Let H be an Ansatz for the expression of h in terms of z[2] and z[3]:
H := c1 * z2^2 + c2 * z3 + c3 * z2 + c4;
IsZero := subs(z2=z[2], z3=z[3], x=t,y=f, h-H):
'h' = eval(H, SolveCoeffsZero(IsZero, q, 40));


# Compute the algebraic relation between z[2] and z[3] using an Ansatz:
H := z3^2 + c1 * z3 + c2 * z2*z3 + c3 * z2^3 + c4 * z2^2 + c5 * z2 + c6:
IsZero := subs(z2=z[2], z3=z[3], x=t, y=f, H):
sort( eval(H,SolveCoeffsZero(IsZero, q, 40)), z3) = 0;
ifactor(j_invariant(lhs(%), z2, z3));