
reset
magma
ZZ:=Integers();
p:=2;
 
F2:=FiniteField(p);

/*
//===== Hermitian curve over F16 ========== // n+1=65 places, genus 6, 3t<n-24=40 =>  t=13 ok => degree(D)=2*6+13=25
r:=4;
F<dt> := FiniteField(p^r); 
q:=p^(ZZ!(r/2));
A<x,T> := AffineSpace(F,2); //doesn't like function fields when asking for an uniformizing parameter at a place
f:=T^q+T-x^(q+1);
K1:=Curve(A,f);
m:=25; //condition for projective normality: m=degree(D) >=2g+1 . Also, t=m-(2g+1)<(n-4g)/3
pole:=true;
Fixed:=[ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 89, 90, 91, 92, 93, 94, 95, 96, 97, 101, 102, 103, 104, 105, 106, 107, 108, 112, 113, 114, 115, 116, 117, 118, 122, 123, 124, 125, 126, 127, 131, 132, 133, 134 ]; // the indices pruned from the _Reduc variables_, not from the _equations_ (because not basis vectors of C2)
*/
/*
//=========== Genus 3 over F32 with record 64 points ============= // ATTENTION: first basis element of L(D) is no more 1 (but fixing it also works, for a reason..)
r:=5;
F<dt> := ext<F2|Polynomial([1,0,1,0,0,1])>;
A<x,y> := AffineSpace(F,2);
X:=x^2+x; Y:=y^2+y;
f:=X^2 + X*Y + dt^3*Y^2 + Y + dt^26;
K1:=Curve(A,f);
m:= 22; //equal 2g+t. with t=17 , such that 3*t<(N-1)-4g=51 -> t< 17 (so no gain to evaluate also at P0..)
pole:=false;
Fixed:=[ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 64, 65, 66, 78 ];
*/

//Elliptic curve over F2^3 with 14 places 
r:=3;
F<dt> := FiniteField(p^r); 
A<x,y> := AffineSpace(F,2);
K1 := Curve(A,y^2+x*y+y-x^3+1);
m:=4; // 2g+t avec t=2 , > 2g+2 = et t=2, verifie bien  6=3t<(n-1)-4g=9
pole:=false;
Fixed:=[];

g:=Genus(K1); g; //
P:=Places(K1,1); #P; //
if pole eq true then n:=#P;
else n:=#P-1; end if;
D:=m*P[1];

Rlist:=[];
for k in [1..101] do
  Append(~Rlist, GaloisRing(p^(k+1), Polynomial([ZZ!Coefficient(MinimalPolynomial(dt,F2),i) : i in [0..r]]) )  );
end for;
R4<Dt>:=Rlist[1];  // $.1^4 + $.1 + 1 for dt  //ResidueField(R4) eq F; //false: doesn't identify the two
R8:=Rlist[2];
ElementToSequence(R4!Dt);



//the evaluation code of L(D) and L(2D) over (P1 if pole) P2,...,P_(n+1)
AGCodes:=function(D,P, pole);
  E:=Basis(D); 
  dimL:=#E; //m+1-Genus
  D2:=2*D;
  E2:=Basis(D2); 
  dimL2:=#E2;
  V,phi:=RiemannRochSpace(D); //(phi^-1)(E[1]); //gives the decomposition of the section E[1] in the basis E of V:  V: (      1       0       0)
  if pole eq true then
    tpole:=UniformizingParameter(P[1]);
    tvp:=tpole^m;
    EvD := Matrix(F,dimL, n, [<i,j,Evaluate(E[i],P[j])>: j in [2..n], i in [1..dimL]]); 
    EvD2 := Matrix(F,dimL2, n, [<i,j,Evaluate(E2[i],P[j])>: j in [2..n], i in [1..dimL2]]);
    for i in [1..dimL] do
      EvD[i,1]:=Evaluate(E[i]*tvp,P[1]);
    end for;
    for i in [1..dimL2] do
      EvD2[i,1]:=Evaluate(E2[i]*tvp^2,P[1]);
    end for;
  else
    EvD := Matrix(F,dimL, n, [<i,j,Evaluate(E[i],P[j+1])>: j in [1..n], i in [1..dimL]]); 
    EvD2 := Matrix(F,dimL2, n, [<i,j,Evaluate(E2[i],P[j+1])>: j in [1..n], i in [1..dimL2]]);
  end if;
  return EvD, EvD2, dimL,dimL2;
end function; 

EvD, EvD2, dimL,dimL2:=AGCodes(D,P,pole);


C:=RowSpace(EvD); //the code of L(D)
dimC:=Dimension(C); //of dimension dimL(D), as expected (here: 3)
CD2:=RowSpace(EvD2); // the code of L(D^2)


//the square of C: dimension checks. 
square:=ZZ!(dimC*(dimC+1)/2);
RowsC2:=[];
for k in [1..dimL] do
  for l in [k..dimL] do
    Append(~RowsC2,  [EvD[k][j]*EvD[l][j] : j in [1..n]]  );
  end for;
end for;

MatC2:=Matrix(F,RowsC2);
C2:=RowSpace(MatC2);
dimC2:=Dimension(C2); 
nbEqVec:=square-dimC2; //square; //square-dimC2; //number of vectorial equations

//extracting a basis of C2 
Bx:=[MatC2[1]]; //
Ix:=[[1,1,1]]; //keeping track of the indexes - and of the ij indexes
C2:=VectorSpaceWithBasis(Bx);
ij:=0;
for i in [1..dimC] do
  for j in [i..dimC] do
    ij+:=1;
    cij:=MatC2[ij];
    if (cij notin C2) then
      Append(~Bx,cij);
      Append(~Ix,[i,j,ij]);
      C2:=VectorSpaceWithBasis(Bx); //f will kill me
    end if;
  end for;
end for;


//decomposition of vectors of C2 on this basis: binary only decompositions! And how sparse !!!
Reduc:=[];
ijsquare:=0;
b:=1;
for i in [1..dimC] do
  for j in [i..dimC] do
    ijsquare+:=1;
    if (i ne Ix[b][1]) or (j ne Ix[b][2]) then
      Append(~Reduc, Coordinates(C2, C2!RowsC2[ijsquare]));
    else
      if b lt dimC2 then 
        b+:=1;
      end if;
    end if;
  end for;
end for;


[i: i in [1..#Reduc] |#[j: j in [1..dimC2] | Reduc[i][j] ne 0] eq 1]; //outputs Fixed


//lift of C
liftGal:=function(a,R); 
  pol:=ElementToSequence(F!a, F2); //
  polZ:=[ZZ!coeff : coeff in pol];
  return R!( &+[polZ[i]*R.1^(i-1): i in [1..r] ] );
end function;
//liftGal(Dt^5,R4);


RowsCR4:=[[liftGal(EvD[i][j],R4) : j in [1..n]]: i in [1..dimC] ] ;


// Square of the lift of C, with basis put aside
CRSquared:=function(CR,Rlift);
  CR2:=[]; BxR:=[]; b:=1;
  for i in [1..dimL] do
    for j in [i..dimL] do
      cij:=  [CR[i][crd]*CR[j][crd] : crd in [1..n]]  ;
      if Ix[b][1] ne i or Ix[b][2] ne j then
        Append(~CR2,cij  );
      else
        Append(~BxR, Vector(cij)  );
        if b lt dimC2 then 
          b+:=1;
        end if;
      end if;
    end for;
  end for;
  MCR2:=Matrix(Rlift,CR2);
  return MCR2,BxR;
end function;
MC24,Bx4:=CRSquared(RowsCR4,R4);


Reduc4:=[Vector([liftGal(Reduc[i][j],R4) : j in [1..dimC2]] ) : i in [1..nbEqVec]];


// Dividing by p^k
dividek:=function(b,R,k);
  pol:=ElementToSequence(R!b);
  polk:=[F2!((ZZ!coeff)/2^(k)) : coeff in pol];
  return F!SequenceToElement(polk,F);
end function;
dividek(4*(R8.1^2+R8.1^3),R8,2); //SequenceToElement([F2|0,0,1,1],F);  // returns dt^6: he converts automatically to the (shorter) primitive element power !


//list of error terms divided by p^k for basis reconstruction of C2 (minus basis) in Z/p^(k+1)Z[X]/P.
errork:=function(MCR2,BxR,ReducR,Rlift,k) 
  pDR:=[];
  for ij in [1..nbEqVec] do
    Append(~pDR,MCR2[ij]-&+[ReducR[ij][l]*BxR[l] : l in [1..dimC2]]);
  end for;
  DF:=[ Vector([dividek(pDR[ij][j],Rlift,k): j in [1..n]]) : ij in [1..nbEqVec]];  //pD4[50]; DF[50]; division seems ok
  return DF;
end function;
DF:=errork(MC24,Bx4,Reduc4,R4,1);


// ======= Multiplicative lift modulo 4 : solve V.A = W in vector V => constructing transp(A)
//transp(W)= DF= transp(A).transp(V) , which expresses on each of the n Big blocks of square-dimC2 rows the vectorial equations: -(ei.fj+ej.fi) + sum_base( mu_(ij,b) (e_bi.e_bj) + Reduc_(ij,b) (e_bi.f_j + e_bj.fi) )
// V: firstly the lifts fi:  n Big blocks of size dimC (lifts of C) for each fixed coord crd. Then, (square-dimC2) Big blocks (coefficients mu_ij for ij fixed) of size dimC2 (coordinates on (e_bi.e_bj)_b for b in IxM)
// DF: as for A, n Big blocks (one per coordinate) of size square

nbVecInc:=dimC-1;// dimC-1;
tA:=ZeroMatrix(F,nbEqVec*n,nbVecInc*n+dimC2*(nbEqVec-#Fixed) ); // (rows: n coords times square vector equations, cols: n blocs of nbVecInc-1 coordinates plus dimC2 coeffs Reduc_ij for each couple ij)
DFvec:=Vector([F!0: j in [1..nbEqVec*n]]);
offset:=nbVecInc*n; // Place taken at the beginning of V by the lifts of C 


populateRow:=procedure(~M,~D,i,j,ijrowA,ijcolA,crd);
//------------POPULATE the invCrd coordinate differently
  row:=(crd-1)*nbEqVec+ijrowA ; //the line of A and D to be changed: inside big block number crd-1 , go fetch couple (i,j) number ij
  idmuij:=offset+dimC2*(ijcolA-1); //skip offset, and fetch the ij-th big block for coefficients mu(ij,l), l varying in dimC2. 
  idfj:=(crd-1)*nbVecInc+j; //go in block crd and look for index of fj
  idfi:=(crd-1)*nbVecInc+i; //go in block crd and look for index of fi
  if j ne 1 then // should remove the 1st column in each of the n first blocks instead
    M[row,idfj]:=-EvD[i][crd];
  end if;
  if i ne 1 then 
    M[row,idfi]:=M[row,idfi]-EvD[j][crd];
  end if;
  D[row]:=DF[ijrowA][crd];
  for b in [1..dimC2] do
    bi,bj:=Explode(Ix[b]); //the indexes i,j of the b-th basis vector ei.ej of C^2
    idbj:=(crd-1)*nbVecInc+bj;
    idbi:=(crd-1)*nbVecInc+bi;
    if bj ne 1 then 
      M[row,idbj]:= M[row,idbj]+ Reduc[ijrowA][b]*EvD[bi][crd] ;
    end if;
    if bi ne 1 then 
      M[row,idbi]:= M[row,idbi]+ Reduc[ijrowA][b]*EvD[bj][crd] ;
    end if;
    if ijrowA notin Fixed then
      M[row,idmuij+b]:= M[row,idmuij+b]+ Bx[b][crd]; //
    end if;
  end for;
end procedure;


ijrowA:=0;
ijcolA:=0;
b:=1;
for i in [1..dimC] do
  for j in [i..dimC] do
    if (i ne Ix[b][1]) or (j ne Ix[b][2]) then
      ijrowA+:=1;
      if ijrowA notin Fixed then
        ijcolA+:=1;
      end if;
      for crd in [1..n] do
        populateRow(~tA,~DFvec,i,j,ijrowA,ijcolA,crd);
      end for;
    else
      if b lt dimC2 then
        b:=b+1;
      end if;
    end if;
  end for;
end for;


//Ker:=NullspaceOfTranspose(tA);
timestart:=Cputime();
Sol,Ker:=Solution(Transpose(tA),DFvec);
print Cputime(timestart);
procedure printInfo()
print "dimC ", dimC, "equals dimL ",dimL, " : ",dimC eq dimL," -- length ",n;
print "dimC2 ", dimC2, " vs total number of unordered couples (i,j): ", square;
print "Mumford's normal generation-- ",dimC2 eq Dimension(CD2);
//print "indexes (k,l) of the extracted basis (ek*el) of C2 ", Ix;
//print "Decomposition of each nonbasis element ei*ej on the previous basis (ek*el) of C2 ", Reduc[[ij : ij in  [1..square] | ij notin [Ix[j][3]: j in [1..dimC2] ] ]];
//print "Error terms, after lifting of this decomposition, for each (i,j)", pD4;
//print "Indexes (i,j) with nonzero error term after lifting of this decomposition ",[i: i in [1..square] |  not IsZero(pD4[i])];
print "Size of the system: #rows of tA ",Nrows(tA)," vs #columns of tA ", Ncols(tA) ;
//print "lignes matrice A non nulles ",[i: i in [1..square*n] |  not IsZero(tA[i])];
//print "indices ij ou A non nulle ",[i: i in [1..square] |  &or[not IsZero(tA[i+k*square]) : k in [0..n-1]] ];
/*print "RowsCR ",RowsCR;
print "RowsC2R ",RowsC2R;
print "RowsCR2 ", RowsCR2; */
//print "==== SOLUTION: ====";
//print "Corrective terms of the lifts of the code + of the lifts of linear relations of C2, that cancel the previous error terms ",Solution(Transpose(tA),DFvec);
print "dimKer ", Dimension(Ker); 
end procedure;
printInfo();


//==========================================================================================
//============= Reconstructing the corrected lift from the solution of the linear system ===
//==========================================================================================

unfoldVects:=function(sol);
  fiVects:=[Vector(F,[sol[(crd-1)*nbVecInc+1] : crd in [1..n] ])];
  for i in [2..dimC] do
    Append(~fiVects,Vector(F,[sol[(crd-1)*nbVecInc+i]: crd in [1..n] ]) );
  end for;
  return fiVects;
end function;


unfoldReduc:=function(sol); //unfolds corrections to the coefficients (=stored in "Reduc") 
  muVects:=[];
  ijcolA:=0;
  for rowA in [1..nbEqVec] do
    if rowA notin Fixed then
      ijcolA+:=1;
      idmuij:=offset+dimC2*(ijcolA-1);
      Append(~muVects,Vector(F,[sol[idmuij+b]: b in [1..dimC2] ]));
    else;
      Append(~muVects,Vector([F!0: b in [1..dimC2]]));
    end if;
  end for;
  return muVects;
end function;


liftGalk:=function(a,R,Rlift); 
  pol:=ElementToSequence(R!a); //! coeff constant a gauche.  Par contre juste Polynomial(dt^3); renvoie renvoie un poly de degre 2, sur F4  ElementToSequence(dt^5, F2);
  polZ:=[ZZ!coeff : coeff in pol];
  return Rlift!( &+[polZ[i]*(Rlift.1)^(i-1): i in [1..r] ] );
end function;
//liftGalk(Dt^5,R4,R8);


liftCorrectC:=function(rowsC,unfoldVects,R,Rlift,k);
  CorrectC:=[Vector(Rlift, [liftGalk(rowsC[1][j],R,Rlift): j in [1..n] ])];
  for i in [2..dimC] do
    Append(~CorrectC,  Vector( [liftGalk(rowsC[i][j],R,Rlift) + p^k*liftGal(unfoldVects[i][j],Rlift  ) : j in [1..n]   ]  )  );
  end for;
  return CorrectC;
end function;


liftCorrectReduc:=function(ReducR,unfoldReduc,R,Rlift,k);
  ReducRlift:=[];
  for row in [1..nbEqVec] do
    Append(~ReducRlift,  Vector( [liftGalk(ReducR[row][b],R,Rlift) + p^k*liftGal(unfoldReduc[row][b],Rlift) : b in [1..dimC2]   ]  )    );
  end for; 
return ReducRlift;
end function;


//Create the total error vector from the list of error vectors
populateD:=procedure(~D,Dlist,ijrowA,crd);
  row:=(crd-1)*nbEqVec+ijrowA ; //the line number of D to be changed: inside big block number crd-1 , go fetch couple (i,j) number ij
  D[row]:=Dlist[ijrowA][crd];
end procedure;


//=====================================================================================
//============= Recursively searching for higher lifts modulo higher p^\ell ===========
//=====================================================================================

findSolk:=function(solR,RowsCR,ReducR,R,Rlift,k);//,time);
  for round in [1] do
    print round;
    solRand:=solR; //+Random(Ker);
    CRlift:=liftCorrectC(RowsCR,unfoldVects(solRand),R,Rlift,k);
    ReducRlift:=liftCorrectReduc(ReducR,unfoldReduc(solRand), R, Rlift,k);
    MCRlift2,BxRlift:=CRSquared(CRlift,Rlift);
    Dk:=errork(MCRlift2,BxRlift,ReducRlift,Rlift,k+1);

    Dkvec:=Vector([F!0: j in [1..nbEqVec*n]]);
    for ijrowA in [1..nbEqVec] do
      for crd in [1..n] do
        populateD(~Dkvec,Dk,ijrowA,crd);
      end for;
    end for;

    if round in [2000*i: i in [1..2^9]] then
      print round;
    end if;

    if IsConsistent(Transpose(tA),Dkvec) then
      solk:=solRand;
      timestart:=Cputime();
      sollift:=Solution(Transpose(tA),Dkvec);
      print Cputime(timestart); //timeend:=Cputime(timestart);
      //timelift:=time+timeend;
      print "found at step ",k;//,sol4,"errorVec ",sol8;
      break;
    end if;
  end for;
  return sollift,CRlift, ReducRlift;//,timelift;
end function;


solR:=Sol; ReducR:=Reduc4; RowsCR:=RowsCR4; R:=R4; Rlift:=Rlist[2]; //time:=0;
for k in [1..99] do
  sollift,CRlift, ReducRlift:=findSolk(solR,RowsCR,ReducR,R,Rlift,k);//,timelift:=    ,time);
  solR:=sollift;
  RowsCR:=CRlift;
  ReducR:=ReducRlift;
  R:=Rlift;
  Rlift:=Rlist[k+2];
  //time:=timelift;
end for;

//=============================================
//============= for the toy example ===========
//=============================================

CorrectC:=[Vector(RowsCR4[1])];
  for i in [2..dimC] do
    Append(~CorrectC,  Vector( [RowsCR4[i][j] + p*liftGal(unfoldVects(Sol)[i][j],R4 ) : j in [1..n]   ]  )  );
  end for;
CorrectC;

CorrectReduc:=[];
  for row in [1..nbEqVec] do
    Append(~CorrectReduc,  Vector( [Reduc4[row][b] + p*liftGal(unfoldReduc(Sol)[row][b],R4) : b in [1..dimC2]   ]  )    );
  end for;

cpp:=function(vec1,vec2);
return Vector([ vec1[i]*vec2[i] : i in [1..n]  ] );
end function;

[&+[  cpp(CorrectC[Ix[b][1]],  CorrectC[Ix[b][2]])*CorrectReduc[i][b] : b in [1..dimC2] ]: i in [1,2]];
cpp(CorrectC[2],  CorrectC[2]); cpp(CorrectC[4],  CorrectC[4]); //The two rows of Reduc (2 = square-dimC2)

