/*

Compute the function j2 where (j0, j1(j0), j2(j0)) is the reduction
modulo p^3 of the j-invariant of the canonical lifting of j0.

*/

load 'witt_gen.m';  
load 'sec_coord.m';
load 'j1.m';

p:=41; // characteristic

P<x0,x1,x2,y0,y1,y2>:=PolynomialRing(GF(p),6);

// to deal with the mod. polyn.
PP<X,Y>:=PolynomialRing(Integers(),2);
PPP<XX,YY>:=PolynomialRing(ResidueClassRing(p^3),2);


R<jj0>:=PolynomialRing(GF(p));
FF<j0>:=FieldOfFractions(R);

// FOR TESTING PURPOSES: compare with what we had computed before by different
// methods.  J1 and J2 were the results we had.
if p eq 5 then
    J1:=3*j0^3 + j0^4;
    J2:=3*j0^5 + 2*j0^10 + 2*j0^13 + 4*j0^14 + 4*j0^15 + 4*j0^16 + j0^17 
        + 4*j0^18 + j0^19 + j0^20 + 3*j0^23 + j0^24;
elif p eq 7 then
    J1:=3*j0^5+5*j0^6;
    J2:=(3*j0^21 + 6*j0^28 + 3*j0^33 + 5*j0^34 + 4*j0^35 + 2*j0^36 
        + 3*j0^37 + 6*j0^38 + 3*j0^39 + 5*j0^40 + 5*j0^41 + 5*j0^42 
        + 2*j0^43 + 3*j0^44 + 6*j0^45 + 3*j0^46 +  5*j0^47 + 5*j0^48 
        + 3*j0^49 + 3*j0^54 + 5*j0^55)/(1 + j0^7);
elif p eq 11 then
    J1:=8*j0^7 + 5*j0^8 + 8*j0^9 + 4*j0^10;
    J2:=(7*j0^33 + 10*j0^44 + 3*j0^55 + 4*j0^66 + 4*j0^73 + 8*j0^74 +
        4*j0^75 + 2*j0^76 + 7*j0^77 + 3*j0^78 + 4*j0^79 + 10*j0^80 + 6*j0^81 +
        10*j0^82 + 8*j0^83 + 3*j0^84 + 2*j0^86 + 4*j0^88 + 7*j0^89 + j0^90 + 
        3*j0^91 + j0^92 + 3*j0^93 + 3*j0^95 + 7*j0^96 + 3*j0^97 + 10*j0^98 + 
        4*j0^99 + 5*j0^100 + 2*j0^101 + 5*j0^102 + 6*j0^103 + 9*j0^104 + 10*j0^105 + 
        7*j0^106 + 4*j0^108 + 9*j0^109 + 7*j0^111 + 4*j0^112 + 4*j0^113 + 9*j0^114 + 
        4*j0^116 + 9*j0^117 + 4*j0^118 + 2*j0^119 + 3*j0^120 + 10*j0^121 + 7*j0^128 + 
        3*j0^129 + 7*j0^130 + 9*j0^131)/(6 + 5*j0^11);
elif p eq 13 then
    J1:=(9*j0^9 + 10*j0^10 + j0^11 + 3*j0^12 + 6*j0^13)/(10 + 11*j0);
    J2:=(10*j0^65 + 11*j0^66 + 2*j0^78 + 10*j0^79 + 12*j0^91 + 8*j0^92 +
        4*j0^104 + 7*j0^105 + 10*j0^113 + j0^114 + 4*j0^115 + 12*j0^116 +
        12*j0^119 + 3*j0^120 + 10*j0^122 + 6*j0^123 + 2*j0^125 + 2*j0^127 + 9*j0^128 + 
        5*j0^129 + 6*j0^130 + 7*j0^131 + 11*j0^132 + 6*j0^133 + 9*j0^134 + j0^135 + 
        12*j0^136 + 7*j0^137 + 2*j0^138 + 12*j0^139 + 11*j0^140 + 2*j0^142 + 9*j0^143 + 
        3*j0^144 + 6*j0^145 + 3*j0^146 + 7*j0^147 + 8*j0^148 + 6*j0^149 + j0^150 + 4*j0^151 + 
        10*j0^152 + 10*j0^153 + 3*j0^154 + 3*j0^155 + 8*j0^156 + 9*j0^157 + 12*j0^158 + 
        7*j0^159 + 3*j0^160 + 7*j0^161 + j0^163 + 7*j0^164 + 9*j0^165 + 9*j0^166 + 10*j0^167 + 
        7*j0^169 + 9*j0^170 + 6*j0^171 + 2*j0^172 + 4*j0^173 + 12*j0^174 + 2*j0^175 + 
        4*j0^176 + 6*j0^177 + 3*j0^178 + 2*j0^179 + 9*j0^180 + 2*j0^181 + 5*j0^182 + 
        4*j0^183 + 10*j0^191 + j0^192 + 4*j0^193 + 12*j0^194 + 11*j0^195)/(12 + 8*j0 + 
    3*j0^13 + 2*j0^14 + j0^26 + 5*j0^27);
elif p eq 17 then
    J1:=(15*j0^11 + 16*j0^12 + 12*j0^13 + 15*j0^14 + 5*j0^15 + 10*j0^16 + 10*j0^17)/(14 + 11*j0);
elif p eq 19 then
    J1:=(j0^13 + 6*j0^14 + 8*j0^15 + 18*j0^16 + 17*j0^17 + 4*j0^18 + 7*j0^19)/(10 + 4*j0);
elif p eq 23 then
    J1:=(19*j0^15 + 22*j0^16 + 12*j0^17 + 13*j0^18 + 19*j0^19 + 15*j0^20 + j0^21 
        + 16*j0^22 + 10*j0^23)/(18 + 16*j0);
elif p eq 29 then
    J1:=(25*j0^19 + 14*j0^20 + 20*j0^21 + 13*j0^22 + 10*j0^23 + 6*j0^24 +
        27*j0^25 + 8*j0^26 + 24*j0^27 + 10*j0^28 + 25*j0^29 + 10*j0^30)/(21 +
    2*j0 + j0^2);
elif p eq 31 then
    J1:=(7*j0^21 + 14*j0^22 + 5*j0^23 + 17*j0^24 + 28*j0^25 + 4*j0^26 +
        12*j0^27 + 7*j0^28 + 5*j0^29 + 8*j0^30 + 17*j0^31)/(27 + 3*j0 +
    15*j0^2);
elif p eq 37 then
    J1:=(29*j0^25 + 5*j0^26 + 28*j0^27 + 33*j0^28 + 32*j0^29 + 9*j0^30 + 9*j0^31 + 
        19*j0^32 + 27*j0^33 + 28*j0^34 + 14*j0^35 + 35*j0^36 + 34*j0^37 + 
        6*j0^38 + 31*j0^39)/(35 + 26*j0 + 16*j0^2 + 20*j0^3);
else
    J1:=FF!fctJ1(p);
end if;



f:=(x0-y0^p)*(x0^p-y0);  // red. mod p of modular pol. (Kronecker rel.)

// partial derivatives (to use the formula for the GT)
fx:=(x0^p-y0);
fy:=(y0^p-x0);
fxx:=0;
fyy:=0;
fxy:=-1;

// modular polynomial (modulo p^3 only is necessary!)
mp:=PP!(PPP!ClassicalModularPolynomial(p));
mpx:=PP!(PPP!Derivative(mp,1,1));
mpy:=PP!(PPP!Derivative(mp,1,2));


// pol = pol in over Z (ONLY WORKS OVER Z!!!)
// returns pol over Z that reduces to sum of a_{i,j,1} x0^{i*p} y0^{j*p}
// Necessary for the formula of the GT
fct1 := function(pol,p)

    P:=Parent(pol);
    ZZ:=Integers();
    v:=Terms(pol);
    res:=0;
    
    for m in v do
        c:=LeadingCoefficient(m);

        if Rank(P) eq 1 then
            nm:=(P.1)^(Degree(m));
        else
            nm:=LeadingMonomial(m);
        end if;
                
        d:=((c - c^p) div p) mod p;
        res+:=d*nm^p;
    end for;

    return res;

end function;


// pol = pol in over Z (ONLY WORKS OVER Z!!!)
// returns pol over Z that reduces to sum of a_{i,j,2} x0^(i*p^2) y0^(j*p^2)
// Needs to write integers as Witt vectors, and so uses 'witt_gen.m'

fct2 := function(pol,p : primr:=[])

    if primr eq []
       then primr:=GenPrimRoots(p,3);
    end if;
    
    P:=Parent(pol);
    ZZ:=Integers();
    v:=Terms(pol);
    res:=0;
    
    for m in v do
        c:=LeadingCoefficient(m);

        if Rank(P) eq 1 then
            nm:=(P.1)^(Degree(m));
        else
            nm:=LeadingMonomial(m);
        end if;

        vc:=IntToWitt(c,p,3 : primr:=primr);
        d:=ZZ!(vc[3]);
        res+:=d*nm^(p^2);
    end for;

    return res;

end function;


/*
pthroot := function(pol)

    P:=Parent(pol);
    p:=Characteristic(P);
    n:=Rank(P);
    v:=Monomials(pol);

    res:=0;
    
    for m in v do
        nm:=1;
        for i in [1..n] do
            deg:=Degree(m,i);
            if ((deg mod p) eq 0) then
                nm*:=(P.i)^(deg div p);
            else
                print "ERROR!!!";
                return -1;
            end if;
        end for;

        res+:=MonomialCoefficient(pol,m)*nm;
    end for;

    return res;

end function;
*/

// This is less generic than the one above but works better here
// Takes p-th root of a polynomial (over Z/(p)) that is a p-th power
pthroot := procedure(~pol)

    P:=Parent(pol);
    p:=Characteristic(P);

    for m in Terms(pol) do
        nm:=1;
        deg:=Degree(m);
        if ((deg mod p) eq 0) then
            nm*:=(P.1)^(deg div p);
        else
            print "ERROR!!!";
            pol:=-1;
            break m;
        end if;

        coef:=LeadingCoefficient(m); // coeff. must be over GF(p)!!!!
        
        pol+:=(coef*nm-m);
    end for;

end procedure;


/*
pthroot := function(pol)

    P:=Parent(pol);
    p:=Characteristic(P);
    v:=Terms(pol);

    res:=0;
    
    for m in v do
        nm:=1;
        deg:=Degree(m);
        if ((deg mod p) eq 0) then
            nm*:=(P.1)^(deg div p);
        else
            print "ERROR!!!";
            return -1;
        end if;
        
        res+:=LeadingCoefficient(m)*nm;
    end for;

    return res;

end function;
*/


// We will leave the term in j2 separated!!
evalf := function(pol)

    return Evaluate(pol,<j0,J1,0,j0^p,J1^p,0>);

end function;

print "Computing pol1x";
pol1x:=Evaluate(fct1(mpx,p),<x0,y0>);

print "Computing pol1y";
pol1y:=Evaluate(fct1(mpy,p),<x0,y0>);

print "Computing pol2";
pol2:=Evaluate(fct2(mp,p),<x0,y0>);

print "Computed pol1x, pol1y, pol2.";

/*
GT:=evalf(fx^(p^2)*x2 + fy^(p^2)*y2 + pol1x^p*x1^p + pol1y^p*y1^p
    + (-1)^(p^2)*x1^p*y1^p + pol2);

GT:=evalf(fy^(p^2)*y2 + pol1x^p*x1^p + pol1y^p*y1^p
    + (-1)^(p^2)*x1^p*y1^p + pol2);
*/

// now without the ``fy^(p^2)*y2'' term

ppol2:=pol1x^p*x1^p + pol1y^p*y1^p
    + (-1)^(p^2)*x1^p*y1^p + pol2;

print "Computed ppol2.";

GT:=evalf(ppol2);

delete pol1x, pol1y, pol2, ppol2, mpx, mpy;

print "Evaluated ppol2";

/*
GT:=evalf(pol1x^p*x1^p + pol1y^p*y1^p
    + (-1)^(p^2)*x1^p*y1^p + pol2);

delete pol1x, pol1y, pol2, mpx, mpy;
*/

// /////////////////////////////////////////////////////

tmp:=Evaluate(fct1(mp,p),<x0,y0>);

delete mp;

/*
g:=fx^p*x1 + fy^p*y1 + tmp;
*/
g:=fy^p*y1 + tmp;
// SpcFct(fx^p*x1 + g) = SpcFct(fx^p*x1) + SpcFct(g) + sum
// The sum is zero after replacing y0 by x0^p in fx, since the summands are multiples of fx.
// SpcFct(fx^p*x1) is also zero (after replacing...); just work with the two terms

print "Computing pol3....";
// pol3:=evalf(SpcFct(g));

// let's do it in a smarter (hopefuly faster) way
EvalSpcFct := function(f)

    v:=Terms(f);
    n:=#v;
    
    if n le 1 then
        return 0;
    end if;

    P:=Parent(f);
    p:=Characteristic(P);

    g1:=&+[ v[i] : i in [1..(n div 2)] ];
    g2:=&+[ v[i] : i in [((n div 2)+1)..n] ];

    gg1:=&+[ evalf(v[i]) : i in [1..(n div 2)] ];
    gg2:=&+[ evalf(v[i]) : i in [((n div 2)+1)..n] ];

    delete v;
    
    tmp:= &+[ -Coef(p,i)*gg1^i*gg2^(p-i) : i in [1..(p-1)] ];

    delete gg1, gg2;
    
    return $$(g1) + $$(g2) + tmp;
    
end function;

pol3:=EvalSpcFct(g);

print "Done!";
    
/*
pol4:=evalf(fSpcFct(g,SpcFct(f))); // should be zero (evalf(SpcFct(f) is zero)
pol5:=evalf(Spc4Fct(f));  // should be zero

if (pol4 ne 0) or (pol5 ne 0) then
    print "check pol4 and pol5";
end if;
*/

// /////////////////////////////////////////////////////

GT+:=pol3;
GT/:=(j0^(p^2)-j0)^(p^2);

delete pol3;
// save "j2_29.magma_session";


print "Finished GT";
print "";
// print GT;
print "";
// ShowMemoryUsage();
// %! free
print "";
print "*************************************************";
print "";

num1:=-Numerator(GT);
den1:=Denominator(GT);

delete GT;

pthroot(~num1);

// ShowMemoryUsage();
// %! free
print "";
print "*************************************************";
print "";

pthroot(~den1);

print "Finished pthroot(Denominator(GT))";
print "";
// ShowMemoryUsage();
// %! free
print "";
print "*************************************************";
print "";


JJ2:=num1/den1;

print "Finished JJ2:=num1/(den1*den)";
print "";
// ShowMemoryUsage();
// %! free
print "";
print "*************************************************";
print "";

print JJ2;


// Test:
if p le 13 then
    JJ2 eq J2;
end if;

