/*

Computes the canonical lifting and elliptic Teichmuller lift of points
of an elliptic curve or the minimal degree lifting.

*/

// NEWEST VERSION!!!!
// Uses the new functions AND uses the whole matrix and let's MAGMA find
// the solution instead of trying to simplify.  -- It will take more memory,
// but maybe not much and seems faster...

// TESTED: works!! -- tested memory and times (compared with old) for p=11 and 13.
//                    (See time_mem_fctWitts.txt.)  Differences seem negligible...


load 'sec_coord.m';

// Computes canonical lift.

    

//  ------- fct ---------

lift := function(a0,b0 : tm:=false, test:=false, ncoord:=2, minimal:=false, choice:=2)
// tm = time it?
// ncoord = no. of coordinates to compute (1 -> x1,P1; 2 -> x2,P2)
// test = test result
// minimal = compute minimal instead of canonical lift?
// choice = make a_1=a_2=0 (choice:=2) or b_1=b_2=0 (choice:=1)
    
if tm then
    ttime:=Cputime();
    ptime:=Cputime();
    print "Testing input.";
end if;

F1:=Parent(a0); // field of definition

error if not F1 cmpeq Parent(b0),
     "ERROR: ",a0, " and ", b0, " are not in the same field";


p:=Characteristic(F1); // characteristic

error if p eq 0,
     "ERROR: Char. of the field is zero!";

error if 4*a0^3+27*b0^2 eq 0,
     "ERROR: singular curve";


if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Create polynomials.";
end if;


if choice eq 1 then
    // we will assume b1=0;
    PR1<x0,c0,c1,a1>:=PolynomialRing(F1,4);
    b1:=0;
else
    // we will assume a1=0;
    PR1<x0,c0,c1,b1>:=PolynomialRing(F1,4);
    a1:=0;
end if;


f:=PR1!(x0^3+a0*x0+b0); // cubic
sf:=SpcFct(f); // sfct (psi_1) of f

// y0^(p-1)
y0pm1:=f^(((p-1) div 2));

// y0^(p+1)
y0pp1:=y0pm1*f;

// Hasse Invariant
HI:=F1!(Coefficient(y0pm1,x0,p-1));

error if HI eq 0,
  "ERROR: hasse inv. = 0!";


if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Compute x1.";
end if;

// dx1/dx0
tmpx1:=PR1!((HI^(-1))*y0pm1-x0^(p-1));

x1:=Integral(tmpx1,x0)+c0+c1*x0^p;

delete tmpx1;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Compute RHS.";
end if;


// RHS - 2nd coord
rhs1:=PR1!((3*x0^(2*p)+a0^p)*x1+a1*x0^p+b1) + sf;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Long division of RHS.";
end if;

// Find canditate to P1 (y1=y0*P1)
// rem1 must be zero!
rhs1*:=F1!(1/2);
P1:=0;
deg1:=Degree(rhs1,x0);
deg2:=(3*((p+1) div 2));

while deg1 ge deg2 do
    lterm:=Coefficient(rhs1,x0,deg1)*x0^(deg1-deg2);
    rhs1-:=lterm*y0pp1;
    P1+:=lterm;
    deg1:=Degree(rhs1,x0);
end while;

rem1:=rhs1;

delete rhs1;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Separate coefficients that must be zero.";
end if;


// separate the coefficients of rem1 -- they must be zero
vcoef:=Coefficients(rem1,x0);

len:=Degree(rem1,x0)+1;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Create matrix and vector.";
end if;


// 3 by len matrix -- 3 variables: (c0,c1,a1) or (c0,c1,b1)
M1:=Matrix(F1,3,len,
     [[F1!(Coefficient(vcoef[j],i+1,1)) : j in [1..len]] : i in [1..3]]);
// print M1;
 
// We now solve the linear system.
VS1:=VectorSpace(F1,len);

// term independent of the variables
// (note: the 4th variable is x0, which does not show up in the coef.)
v1:=VS1![-Evaluate(vcoef[i],<0,0,0,0>) : i in [1..len]];


if tm 
   then print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Solve the linear system.";
end if;


//atime:=Cputime();
// find ONE solution!
vsol:=Solution(M1,v1);
// print "Time: ", Cputime(atime);

// IMPROVE?
// We could find all...  Check Nullspace
// file:///usr/local/magma/doc/html/text616.htm#5579


// Now, compute the solutions

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Find x1 from solution.";
end if;


// x1 only involves c0 and c1, but not b1
x1:=Evaluate(x1,<x0,vsol[1],vsol[2],0>);

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Find P1 from solution.";
end if;

P1:=Evaluate(P1,<x0,vsol[1],vsol[2],vsol[3]>);

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
end if;

// ********** test ***************

if test then 
    print "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%";
    if Evaluate(rem1,<x0,vsol[1],vsol[2],vsol[3]>) eq 0
	then print "test1 is OK";
    else print "test 1 in bad!"; print "rem1 = ", Evaluate(rem1,<x0,vsol[1],vsol[2],vsol[3]>);
    end if;
    print "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%";
end if;


delete M1, rem1, VS1, vcoef, len, v1;


if tm then 
    print "Done with 1st coord.";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
end if;

// return a1, b1, x1, P1
// if ncoord=1, stops here!
if ncoord eq 1 then
    if choice eq 1 then
        return vsol[3], 0, x1, P1;
    else
        return 0, vsol[3], x1, P1;
    end if;
end if;

//////////////////////////////////////////////////
////////////// SECOND COORDINATE /////////////////
//////////////////////////////////////////////////

//  Now to red. mod p^3

if tm then
    ptime:=Cputime();
    print "Creating new structures.";
end if;

// introduce the new variables: a2 (and b2=0) or b2 (and a2=0)
//  and the coeff. of x0^(pk) in x2
// NOTE: need one less var. than if a2!=0 and b2!=0!!!
if choice eq 1 then
    PR2<x0,a2>:=PolynomialRing(F1,((3*p+5) div 2));
    b1:=0;
    b2:=0;
else
    PR2<x0,b2>:=PolynomialRing(F1,((3*p+5) div 2));
    a1:=0;
    a2:=0;
end if;
  
// We need to convert the previous results to this new ring.
Convert := procedure(~P,M) // convert P to new magma M
    P:= &+[F1!(Coefficient(P,1,i))*(M.1)^i : i in [0..Degree(P,1)]];
end procedure;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Converte polynomials.";
end if;


// make the conversions
Convert(~P1,PR2);
Convert(~x1,PR2);
Convert(~y0pm1,PR2);
Convert(~y0pp1,PR2);
Convert(~f,PR2);
Convert(~sf,PR2);

// set the values of coefficients
if choice eq 1 then
    va1:=vsol[3]; // a1
    vb1:=0; // b1
else
    vb1:=vsol[3]; // b1
    va1:=0; // a1
end if;

delete vsol;

if tm 
   then print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Compute new polynomials.";
end if;


// compute new powers of y0:
// y0^(p^2-1) and y0^(p^2+1)
y0psm1:=y0pm1^p*y0pm1;
y0psp1:=f*y0psm1;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Compute x2.";
end if;

// dx2/dx0
tmpx2:=PR2!((HI^(-1-p))*y0psm1-x0^(p^2-1)-x1^(p-1)*(HI^(-1)*y0pm1-x0^(p-1)));

// we ARE adding the term in x0^(p^2)!!
x2:=Integral(tmpx2,x0)+ &+[PR2.(i+3)*x0^(i*p) : i in [0..((3*p-1) div 2)] ];


delete tmpx2, HI;

// add extra terms (known)
if not minimal then
    tmp:= F1!(3/4)*x1^2;
    x2+:= &+[(Coefficient(tmp,x0,(i div p+p)))^p*x0^i : i in [((3*p^2+p) div 2) .. 2*p^2-p by p]];
    delete tmp;
end if;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Create vector for RHS.";
end if;

// find the necessary vector for the call of SpcFctv
if choice eq 1 then // b1=0
    if a0 eq 0 then
        ev:=[ 3*x0^(2*p)*x1 , va1*x0^p ];
    else
        ev:=[ 3*x0^(2*p)*x1 , a0^p*x1, va1*x0^p ];
    end if;
else // a1=0
    if a0 eq 0 then
        ev:=[ 3*x0^(2*p)*x1 , vb1 ];
    else
        ev:=[ 3*x0^(2*p)*x1 , a0^p*x1, vb1 ];
    end if;
end if;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "First part of RHS.";
end if;


// RHS
rhs2:= (3*x0^(2*p^2)+a0^(p^2))*x2 + 3*x0^(p^2)*x1^(2*p) +
    (-NextCoord(3,p)*x0^(2*p^2)+va1^p)*x1^p + a2*x0^(p^2) + b2;


if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Second part (SpcFct(ev)).";
end if;


rhs2+:= SpcFctv(ev);

delete ev;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Third part (fSpcFct).";
end if;


rhs2+:= fSpcFct( (3*x0^(2*p)+a0^p)*x1 + va1*x0^p + vb1 , sf );

delete sf;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Fourth part (Spc4Fct).";
end if;

rhs2+:= Spc4Fct(f);
    
if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Subtract known terms.";
end if;

// now subtract parts of the LHS that are known!
// what is left is 2*y0^(p^2+1)*P2!
rhs2-:=(f^p*P1^(2*p)+F1!((2-2^p)/p)*y0pp1^p*P1^p); // NOT really rhs2...

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Long division of RHS.";
end if;

// rem2 has to be zero
rhs2*:=F1!(1/2);
P2:=0;
deg1:=Degree(rhs2,x0);
deg2:=(3*((p^2+1) div 2));

while deg1 ge deg2 do
    lterm:=Coefficient(rhs2,x0,deg1)*x0^(deg1-deg2);
    rhs2-:=lterm*y0psp1;
    P2+:=lterm;
    deg1:=Degree(rhs2,x0);
end while;

rem2:=rhs2;

delete rhs2, y0pm1, y0pp1, y0psm1, y0psp1;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Create vector for RHS.";
end if;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Collect terms that must be zero.";
end if;

vcoef:=Coefficients(rem2,x0);

len:=Degree(rem2,x0)+1;


if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Compute matrix and vector.";
end if;

r:=((3*p+3) div 2);  // number of variables (except x0) = rank!

// Use echelon form to find solution
tzero:=<0 : i in [1..(3*p+5) div 2]>;


// CAREFUL HERE!!!
// MAGMA solves system like x*A=b, not A*x=b

// full matrix and vector
M2:=Matrix(F1,r,len,[[F1!(Coefficient(vcoef[j],i+1,1)) : j in [1..len]] : i in [1..r]]);
v2:=Vector([F1!(-Evaluate(vcoef[j],tzero)) : j in [1..len]]);

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Solve the system.";
end if;


vsol:=Solution(M2,v2);


// IMPROVE?
// Find all solutions??


delete M2, v2;

// no. of varibales - 1 (for x0)
// r:=((3*p+3) div 2);

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Create vector to evaluate solutions.";
end if;

// we now create a vector to evaluate the solutions
// it must have <x0,a2 or b2, ..coefs of x2.. >
tp:=<x0>;

for i in [1..r] do
    tp:=Append(tp,vsol[i]);
end for;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
end if;


// ********** test ***************
if test then
    print "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%";
    if Evaluate(rem2,tp) eq 0
	then print "test2 is OK";
    else print "test 2 in bad!"; print "rem2 = ",  Evaluate(rem2,tp);
    end if;
    print "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%";
end if;

delete rem2;

if tm then 
    ptime:=Cputime();
    print "Find x2 from solution.";
end if;

x2:=Evaluate(x2,tp);

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
    ptime:=Cputime();
    print "Find P2 from solution.";
end if;

P2:=Evaluate(P2,tp);

delete tp, r;

if tm then 
    print "Done!";
    print "Partial time = ", Cputime(ptime);
    print "Total time = ", Cputime(ttime);
    print "**************************************************";
end if;


if choice eq 1 then
    return va1, vb1, x1, P1, vsol[1], 0, x2, P2;
else
    return va1, vb1, x1, P1, 0, vsol[1], x2, P2;
end if;


end function;



/* TESTED!!!

boole[~/tex/comp]$ magma
Magma V2.15-4     Fri Jan 23 2009 11:45:37 on boole    [Seed = 2744327905]
Type ? for help.  Type <Ctrl>-D to quit.
> load 'fctWitt2.m';
Loading "fctWitt2.m"
Loading "sec_coord.m"
> load 'fctWitt.m';
Loading "fctWitt.m"
Loading "sec_coord.m"
> p:=5; F<a0,b0>:=RationalFunctionField(GF(p),2); time a1,b1,x1,P1,a2,b2,x2,P2:=lift(a0,b0);
Time: 0.380
> aa1,bb1,xx1,PP1,aa2,bb2,xx2,PP2:=lift2(a0,b0);
> time aa1,bb1,xx1,PP1,aa2,bb2,xx2,PP2:=lift2(a0,b0);
Time: 0.250
> [a1,b1,x1,P1,a2,b2,x2,P2] eq [aa1,bb1,Parent(x1)!xx1,Parent(P1)!PP1,aa2,bb2,Parent(x2)!xx2,Par\
ent(P2)!PP2];
true
> p:=7; F<a0,b0>:=RationalFunctionField(GF(p),2); time a1,b1,x1,P1,a2,b2,x2,P2:=lift(a0,b0);
Time: 2.940
> time aa1,bb1,xx1,PP1,aa2,bb2,xx2,PP2:=lift2(a0,b0);
Time: 2.340
> [a1,b1,x1,P1,a2,b2,x2,P2] eq [aa1,bb1,Parent(x1)!xx1,Parent(P1)!PP1,aa2,bb2,Parent(x2)!xx2,Par\
ent(P2)!PP2];
true
> p:=11; F<a0,b0>:=RationalFunctionField(GF(p),2); time a1,b1,x1,P1,a2,b2,x2,P2:=lift(a0,b0);
Time: 112.630
> time aa1,bb1,xx1,PP1,aa2,bb2,xx2,PP2:=lift2(a0,b0);
Time: 96.550
> [a1,b1,x1,P1,a2,b2,x2,P2] eq [aa1,bb1,Parent(x1)!xx1,Parent(P1)!PP1,aa2,bb2,Parent(x2)!xx2,Par\
ent(P2)!PP2];
true

*/



/*

Test it!!!  Compared the result of this file with some old ones...

p:=7;F<a0,b0>:=RationalFunctionField(GF(p),2);


load 'fctWitt4.m';
a1, b1, x1, P1, a2, b2, x2, P2 := lift(a0,b0 : tm:=true, test:=true, ncoord:=2, minimal:=false, choice:=1);


load 'fctWitt3.m';
aa1, bb1, xx1, PP1, aa2, bb2, xx2, PP2 := lift(a0,b0 : tm:=true, test:=true, ncoord:=2, minimal:=false);

Px1:=Parent(x1);
Px2:=Parent(x2);
a1 eq aa1;
b1 eq bb1;
x1 eq Px1!xx1;
P1 eq Px1!PP1;
a2 eq aa2;
b2 eq bb2;
x2 eq Px2!xx2;
P2 eq Px2!PP2;

load 'fctWitt4.m';
a1, b1, x1, P1, a2, b2, x2, P2 := lift(a0,b0 : tm:=true, test:=true, ncoord:=2, minimal:=false, choice:=2);


load 'fctWitt2.m';
aa1, bb1, xx1, PP1, aa2, bb2, xx2, PP2 := lift(a0,b0 : tm:=true, test:=true, ncoord:=2, minimal:=false);

Px1:=Parent(x1);
Px2:=Parent(x2);
a1 eq aa1;
b1 eq bb1;
x1 eq Px1!xx1;
P1 eq Px1!PP1;
a2 eq aa2;
b2 eq bb2;
x2 eq Px2!xx2;
P2 eq Px2!PP2;



*/







/*

// Let's use a separate function to compute the RHS of the 3rd coord
F3rdCoord := function(v,p)
    // v is a vector which can give us some zero coeff.
    // v=[a0,a1,a2,b0,b1,b2] : a value equal to zero makes the fct
    // assume the the corresponding term is zero.
    F<A0,B0>:=RationalFunctionField(GF(p),2);
    PR<A1,A2,B1,B2,X0,X1,X2>:=PolynomialRing(F,7);

    if v[1] eq 0 then
        a0:=0;
    else
        a0:=A0;
    end if;

    if v[2] eq 0 then
        a1:=0;
    else
        a1:=A1;
    end if;

    if v[3] eq 0 then
        a2:=0;
    else
        a2:=A2;
    end if;

    if v[4] eq 0 then
        b0:=0;
    else
        b0:=B0;
    end if;

    if v[5] eq 0 then
        b1:=0;
    else
        b1:=B1;
    end if;

    if v[6] eq 0 then
        b2:=0;
    else
        b2:=B2;
    end if;

    term1:= (3*X0^(p^2)+a0^(p^2))*X2 + 3*X0^(p^2)*X1^(2*p) +
        (-NextCoord(3,p)*X0^(2*p^2)+a1^p)*X1^p + a2*X0^(p^2) + b2;

    term2:= SpcFct( (3*X0^p+a0^p)*X1 + a1*X0^p + b1 );

    term3:= fSpcFct( (3*X0^p+a0^p)*X1 + a1*X0^p + b1 , SpcFct(X0^3+a0*X0+b0) );

    term4:= Spc2Fct(X0^3+a0*X0+b0);

    term5:= Spc3Fct(X0^3+a0*X0+b0);

    return term1 + term2 + term3 + term4 + term5;

end function;

*/

