/*

Checks j1 and j2 with

  mpol_p((j0,j1,j2),(j0^p,j1^p,j2^p)) = 0

*/

load 'witt_gen.m';
load 'witt2_fun.m';

p:=17;

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];
// vy:=[j0^p,j1^p,j2^p,j3^p];

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

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

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

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

// don't use EvalPolWitt as we used the optimized add and prod.
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:=WittProd2(tmp,vxs[degX+1]);
    tmp:=WittProd2(tmp,vys[degY+1]);
    
    res:=WittSum2(res,tmp);
    
end for;

res; // should be [0,0,0]
