/*

Computes j3 using the modular polynomial (only).

*/

load 'witt_gen.m';

p:=11;

vecS:=WittSumPol(3,p);
vecP:=WittProdPol(3,p);
primr:=GenPrimRoots(p,4);

F<j0>:=RationalFunctionField(GF(p));
P<j3>:=PolynomialRing(F);
PP<X,Y>:=PolynomialRing(Integers(),2);

load 'load_js.m';

mp:=PP!ClassicalModularPolynomial(p);


vx:=[j0,j1,j2,j3];
// vy:=[j0^p,j1^p,j2^p,j3^p];

vxs:=<[F!1,F!0,F!0,F!0], vx>;

// compute powers of vx
for i in [2..p+1] do
    Append(~vxs,WittProd(vxs[i],vx : vecP:=vecP));
end for;

// compute powers of vy
vys:=< WittFrob(v) : v in vxs >;

res:=[F!0,F!0,F!0,F!0];

// evaluate mod pol at ([j0,j1,j2,j3], [j0^p,j1^p,j2^p,j3^p])
for t in Terms(mp) do
    coef:=LeadingCoefficient(t);
    degX:=Degree(t,1);
    degY:=Degree(t,2);
    tmp:=IntToWitt(coef,p,4 : primr:=primr);
    tmp:=WittProd(tmp,vxs[degX+1] : vecP:=vecP);
    tmp:=WittProd(tmp,vys[degY+1] : vecP:=vecP);
    
    res:=WittSum(res,tmp : vecS:=vecS);
    
end for;

eqt:=res[4];
delete res;

den:=Coefficient(eqt,p);
num:=-Coefficient(eqt,0);

j3:=InvFrob(num)/InvFrob(den);
