/*

Finds degrees of the "partial derivatives terms" showing up in the
formula for J2

Also finds the valuations at p of the a0, a1, where the Modular
Polynomial is: a0 + a1(X + Y) + ...

*/

load 'witt_gen.m';

// 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;




for p in [ i : i in [5..59] | IsPrime(i) ] do

    P<X,Y>:=PolynomialRing(Integers(),2);
    mp:=P!ClassicalModularPolynomial(p);

    fx:=Evaluate(fct1(Derivative(mp,1),p),<X,X^p>);
    fy:=Evaluate(fct1(Derivative(mp,2),p),<X,X^p>);

    g:=Evaluate(fct2(mp,p),<X,X^p>);

    v1:=[ Degree(fx,1), Degree(fy,1), Degree(g,1) ];
    v2:=[ p^3+p^2-p, p^3, p^4+p^3-p^2 ];

    // print <p, v1,v2  >; 
    print <p, [ v1[i]-v2[i] : i in [1..3] ]  >; 

    PP<x,y>:=PolynomialRing(GF(p),2);
    
    v1:=[ Degree(PP!fx,1), Degree(PP!fy,1), Degree(PP!g,1) ];
    v2:=[ p^3+p^2-p, p^3-p^2, p^4+p^3-p^2 ];

    // print <p, v1,v2  >; 
    print <p, [ v1[i]-v2[i] : i in [1..3]]  >; 

end for;

/*

Result:

<5, [ 0, 0, 0 ]>
<5, [ 0, 0, 0 ]>
<7, [ 0, 0, 0 ]>
<7, [ 0, 0, 0 ]>
<11, [ 0, 0, 0 ]>
<11, [ 0, 0, 0 ]>
<13, [ 0, 0, 0 ]>
<13, [ 0, 0, 0 ]>
<17, [ 0, 0, 0 ]>
<17, [ 0, 0, 0 ]>
<19, [ 0, 0, 0 ]>
<19, [ 0, 0, 0 ]>
<23, [ 0, 0, 0 ]>
<23, [ 0, 0, 0 ]>
<29, [ 0, 0, 0 ]>
<29, [ 0, 0, 0 ]>
<31, [ 0, 0, 0 ]>
<31, [ 0, -31, 0 ]>
<37, [ 0, 0, 0 ]>
<37, [ 0, 0, 0 ]>
<41, [ 0, 0, -1681 ]>
<41, [ 0, 0, -1681 ]>
<43, [ 0, 0, 0 ]>
<43, [ 0, 0, 0 ]>
<47, [ 0, 0, 0 ]>
<47, [ 0, 0, 0 ]>
<53, [ 0, 0, 0 ]>
<53, [ 0, 0, 0 ]>
<59, [ 0, 0, 0 ]>
<59, [ 0, 0, 0 ]>


> Factorization(1681);
[ <41, 2> ]


*/

/*
myFac := function(n)
    if n eq 0 then
        return <0>;
    else
        return Factorization(n);
    end if;

end function;
*/


// Now, let's see what are the constant terms
for p in [ i : i in [2..59] | IsPrime(i) ] do

    P<X,Y>:=PolynomialRing(Integers(),2);
    mp:=P!ClassicalModularPolynomial(p);
    mpx:=Derivative(mp,1);

    // print <p, myFac(Evaluate(mp,<0,0>)), myFac(Evaluate(mpx,<0,0>)) >;
    <p,  Valuation(Evaluate(mp,<0,0>),p), Valuation(Evaluate(mpx,<0,0>),p), p mod 6 >;
    
end for;

/*

RESULT:

<2, 12, 8, 2>
<3, Infinity, 3, 3>
<5, 3, 3, 5>
<7, Infinity, Infinity, 1>
<11, 3, 3, 5>
<13, Infinity, Infinity, 1>
<17, 3, 3, 5>
<19, Infinity, Infinity, 1>
<23, 3, 4, 5>
<29, 3, 3, 5>
<31, Infinity, Infinity, 1>
<37, Infinity, Infinity, 1>
<41, 3, 3, 5>
<43, Infinity, Infinity, 1>
<47, 3, 3, 5>
<53, 3, 3, 5>
<59, 3, 3, 5>


*/



/*

Ran this by hand: files phi_<prime>.2

> p:=61; load 'phi_61.2';
Loading "phi_61.2"
> <p, Valuation(a0,p), Valuation(a1,p), p mod 6>;

RESULT

<61, Infinity, Infinity, 1>
<67, Infinity, Infinity, 1>
<71, 3, 3, 5>
<73, Infinity, Infinity, 1>
<79, Infinity, Infinity, 1>
<83, 3, 3, 5>
<89, 3, 3, 5>
<97, Infinity, Infinity, 1>
<101, 3, 3, 5>
<103, Infinity, Infinity, 1>
<107, 3, 3, 5>
<109, Infinity, Infinity, 1>
<113, 3, 3, 5>
<127, Infinity, Infinity, 1>
<131, 3, 3, 5>
<137, 3, 3, 5>
<139, Infinity, Infinity, 1>
<139, Infinity, Infinity, 1>
<149, 3, 3, 5>
<151, Infinity, Infinity, 1>
<157, Infinity, Infinity, 1>
<163, Infinity, Infinity, 1>
<167, 3, 3, 5>
<173, 3, 3, 5>
<179, 3, 3, 5>
<181, Infinity, Infinity, 1>
<191, 3, 3, 5>
<193, Infinity, Infinity, 1>
<197, 3, 3, 5>
<199, Infinity, Infinity, 1>




*/


/*

// DOES NOT WORK!!!!!
// for p>59, use file from Rubinstein phi_<prime>.  The phi_<prime>.2
// contains only the coefficiens that we care.
for p in [ i : i in [61..200] | IsPrime(i) ] do

    file:="phi_" cat IntegerToString(p) cat ".2";
    load file;
    
    <p,  Valuation(a0,p), Valuation(a1,p) >;
    
end for;
*/
