/*

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

*/

load 'fctWitt.m';
load 'wittj.m';

fctJ1 := function(p)

    P<j0>:=RationalFunctionField(GF(p));
    
    // precomputed:
    if p eq 5 then
        return 3*j0^3 + j0^4;
    elif p eq 7 then
        return 3*j0^5+5*j0^6;
    elif p eq 11 then
        return 8*j0^7 + 5*j0^8 + 8*j0^9 + 4*j0^10;
    elif p eq 13 then
        return (9*j0^9 + 10*j0^10 + j0^11 + 3*j0^12 + 6*j0^13)/(10 + 11*j0);
    elif p eq 17 then
        return (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
        return (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
        return (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
        return (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
        return (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
        return (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);
    end if;


    if p lt 61 then // precomp. Mod. Pol. in Magma
        PP:=PolynomialRing(ResidueClassRing(p^2),2);
        PPP:=PolynomialRing(Integers(),2);
        MP:=Evaluate(PPP!(PP!ClassicalModularPolynomial(p)),<PPP.1,PPP.1^p>) div p;

        return -Evaluate(MP,<j0,0>)/(j0^(p^2)-j0);

    else // no  precomp. Mod. Pol. in Magma
        PP<a0,b0>:=RationalFunctionField(GF(p),2);
        a1, b1:= lift(a0,b0 : ncoord:=1);
        a:=[a0,a1]; b:=[b0,b1];

        return wittj(a,b : jout:=true);

    end if;
    
        
end function;
