/*

Computes j-inv as Witt vector from coef (vectors) a and b

*/

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

wittj := function(a,b : jout:=false)

    // jout:=true is for the case when the coefficients are given in
    //    F_p(a0,b0) to return the vector in terms of j0 instead of a0,b0
    
    F:=Parent(a[1]);
    p:=Characteristic(F);
    l:=#a;

    primr:=GenPrimRoots(p,l);
    v4:=IntToWitt(4,p,l : primr:=primr);
    v27:=IntToWitt(27,p,l : primr:=primr);
    v1728:=IntToWitt(1728,p,l : primr:=primr);

    
    if l eq 3 then

        v4a3:=WittProd2(v4,WittProd2(a,WittProd2(a,a)));
        v27b2:=WittProd2(v27,WittProd2(b,b));

        num:=WittProd2(v1728,v4a3);
        den:=WittInv(WittSum2(v4a3,v27b2));

        j:=WittProd2(num,den);

    elif l eq 2 then

        PP<X0,X1,Y0,Y1>:=PolynomialRing(Integers(),4);
        
        vecP:=[X0*Y0, X0^p*Y1 + Y0^p*X1];
        vecS:=[X0 + Y0, X1+Y1 - &+[ (Binomial(p,i) div p)*X0^i*Y0^(p-i) :
            i in [1..(p-1)] ] ];

        v4a3:=WittProd(v4,WittPwr(a,3 : vecP:=vecP) : vecP:=vecP);
        v27b2:=WittProd(v27,WittPwr(b,2 : vecP:=vecP) : vecP:=vecP);

        num:=WittProd(v1728,v4a3 : vecP:=vecP);
        den:=WittInv(WittSum(v4a3,v27b2 : vecS:=vecS));

        j:=WittProd(num,den : vecP:=vecP);
        

    else

        vecP:=WittProdPol(l-1,p);
        vecS:=WittSumPol(l-1,p);

        v4a3:=WittProd(v4,WittPwr(a,3 : vecP:=vecP) : vecP:=vecP);
        v27b2:=WittProd(v27,WittPwr(b,2 : vecP:=vecP) : vecP:=vecP);

        num:=WittProd(v1728,v4a3 : vecP:=vecP);
        den:=WittInv(WittSum(v4a3,v27b2 : vecS:=vecS));

        j:=WittProd(num,den : vecP:=vecP);
        

    end if;

    if jout then
        // check if the field is F_p(a0,b0)
        if (Rank(F) eq 2) and (BaseField(F) eq GF(p)) then

            FF<j0>:=RationalFunctionField(GF(p));

            res:=[j0];

            for i in [2..l] do

                ji:=j[i];
            
                num:=0;
                den:=0;

                // we first use j0=a0^3/b0^2
                for term in Terms(Numerator(ji)) do
                    coef:=LeadingCoefficient(term);
                    deg:=Degree(term,1) div 3;
                    num+:=coef*j0^deg;
                end for;

                for term in Terms(Denominator(ji)) do
                    coef:=LeadingCoefficient(term);
                    deg:=Degree(term,1) div 3;
                    den+:=coef*j0^deg;
                end for;

                // now convert to proper: a0^3/b0^2 = 27*j0/(4*(1728-j0))
                num:=Evaluate(num, 27*j0/(4*(1728-j0)));
                den:=Evaluate(den, 27*j0/(4*(1728-j0)));

                Append(~res,num/den);

            end for;

            return res;
            
        else // wrong field

            print "Need field to be GF(p)(a0,b0).";
            print "Field used: ", F;
            
        end if;

    else // jout eq false

        return j;

    end if;
            

end function;

    
