//***************************************************
//Cariter-Manin_CheckRadical.txt
//For rational prime p, 
//Compute Cartier-Manin matrix M of DCEC, and Let I be the ideal generated by entries of M and f2,
//where f2 is given in main theorem1. 
//Check whether I is the radical by using the built-in function "IsRadical".  
//***************************************************

//System("date");
p:=11;
n:=2;
Kprime:=GF(p);
K:=GF(p^n);
INT:=Integers();

CR<a0,a1,a2,a3,a4,a5>:=PolynomialRing(K,6,"grevlex");
R<X,Y,A,B>:=PolynomialRing(CR,4,"lex");
CR_T<t>:=PolynomialRing(K);
PR<X,Y,Z,W>:=PolynomialRing(CR,4,"grevlex");
CR2<a0,a1,a2,a3,a4,a5,s>:=PolynomialRing(K,7,"grevlex");

 
//compute polynomial f2 given in Main Theorem 3.1
NonsingPolynomial:=function();
  CRprime<a0,a1,a2,a3,a4,a5>:=PolynomialRing(Kprime,6,"grevlex");
  Rprime<x,y,A,B>:=PolynomialRing(CRprime,4,"grevlex");
  Sprime<X,Y,Z,AA,BB>:=PolynomialRing(CRprime,5);
  PP := Y^2*Z-(X^3+AA*X*Z^2+BB*Z^3);
  qq := a0*X^2 + a1*X*Y + a2*X*Z^2 + a3*Y^2*Z + a4*Y*Z + a5*Z^2;
  phi := hom<Sprime->Rprime | x,y,1,A,B>;
  P := phi(PP);
  q := phi(qq);
  Res := Resultant(P,q,y);
  nonsing:=[Rprime!0,Rprime!0,Rprime!0];
  
 Start_D:=Realtime();
  DiscG := Discriminant(Res,x);
  FactDiscG := Factorization(DiscG);
  nonsing[1]:= 4*A^3+27*B^2;
  nonsing[3]:= FactDiscG[1][1];
  nonsing[2]:= DiscG div (nonsing[3] ^ (FactDiscG[1][2]));
 End_D:=Realtime();
  printf"Compute time Discriminant= %o\n",End_D - Start_D;

  CR<a0,a1,a2,a3,a4,a5>:=PolynomialRing(K,6,"grevlex");
  R<X,Y,A,B>:=PolynomialRing(CR,4,"lex");
  NonSing:=[];
  for i in [1..3] do 
    NonSing[i] := R!nonsing[i];
  end for;

  return NonSing;
end function;


//for charactaristic p,compute coefficients of supersingular curves
SUPERSINGULAR:=function(p);
  CR_T<t>:=PolynomialRing(K);
  SP:=SupersingularPolynomial(p);
  j_list := RootsInSplittingField(SP);
  model_list := [WeierstrassModel(EllipticCurveFromjInvariant(K!j[1])) : j in j_list];
  coef_list := [Coefficients(w) : w in model_list];
  return coef_list;
end function;

//for charactaristic p,compute coefficients of supersingular curves by Algorithm5.1
CartierManinMtrx:=function(c);
  CR<a0,a1,a2,a3,a4,a5>:=PolynomialRing(K,6,"grevlex");
  PR<X,Y,Z,W>:=PolynomialRing(CR,4,"grevlex");
  P := Y^2*Z - (X^3 + c[4]*X*Z^2 + c[5]*Z^3); // Elliptic curve;  
  //"Elliptic Cueve",P;//
  q:=a0*X^2 + a1*X*Y + a2*X*Z + a3*Y^2 +a4*Y*Z + a5*Z^2;
  Q :=  W^2 -q; //double cover of P   
  F1:=P^(p-1);
  C:=ZeroMatrix(CR,4,4);
  E:=[[1,1,1,2],[1,1,2,1],[1,2,1,1],[2,1,1,1]];
  index:=[0,0,0,0]; //index of x,y,z,w in t

  S_Alltime:=0;
  for i in [1,2,3,4] do
     for j in [1,2,3,4] do
       for k in [1,2,3,4] do
         index[k]:=p*E[i][k] - E[j][k];
       end for;//k   
       if IsDivisibleBy(index[4],2) eq true then
          Start_S:=Realtime();   
           S:=Binomial(p-1, INT!(index[4]/2))*W^index[4]*(-q)^(p-1-INT!(index[4]/2));
           F:=F1*S;
          End_S:=Realtime();
          S_Alltime:=  S_Alltime + End_S - Start_S;
       else 
          F:=PR!0; //S=0, and  F:=P^(p-1)*S;
       end if;
       
       C[i,j] := MonomialCoefficient(F,X^index[1]*Y^index[2]*Z^index[3]*W^index[4]);
    end for;//j
  end for;//i
  printf"compute time Cartier-Manin matrix =%o\n",S_Alltime;
  return C;
end function;

//check (1,1)-entry of Cartier-Manin Matrix is zero (cf. Lemma2.13)
ZeroJudge:=procedure(c);
  if   CR!c ne CR!0 then
    printf"note:Carter-Manin[1,1] = %o\n",c;
  end if;
end procedure;

//change the matrix to the sequence
MtxToSeq:=function(C);
  rho1:=hom<PR -> CR | 1,1,1,1>;
  B:=[];
  n:=INT!1;
  for i in [1,2,3,4] do
     for j in [1,2,3,4] do
      if PR!C[i,j] ne PR!0 then
        B[n]:=rho1(C[i,j]);
        n := n+1;
      end if;
     end for;
  end for;
  return B;
end function;



//-----------------MAIN--------------------------
Start:=Realtime();
printf " ===========K = GF(%o^%o)==============\n",p,n; 
//"charcteristic of base field K:",p;
NS:=NonsingPolynomial();
List_ssingEll:=SUPERSINGULAR(p);

Vtime:=0;
L0:=K;
Varieties:=[];
for i in [1..#List_ssingEll] do 
  coeff_P := List_ssingEll[i];
  A1:=coeff_P[4];
  B1:=coeff_P[5];
  Varieties[i]:=[];
  printf" ==== supsersing. elliptic curve  P := Y^2*Z - (X^3 + %o*X*Z^2 + %o*Z^3)====\n",A1,B1;
  M:=CartierManinMtrx(coeff_P);
  ZeroJudge(M[1,1]); //->nesessary? 
  Basis:=MtxToSeq(M);  
  //-------------- Map from CR to CR2--------------------
  rho2:=hom<R -> CR | 1,1,A1,B1>;
  rho3:=hom<CR -> CR2 | a0,a1,a2,a3,a4,a5>; 
  NS1:=rho2(NS[1]);
  NS2:=rho2(NS[2]);
  NS1:=rho3(NS1);
  NS2:=rho3(NS2); 
  NewBasis:=[];
  for i in [1..#Basis] do
    NewBasis[i] := rho3(Basis[i]);
  end for;
  NewBasis:=Append(NewBasis,s*NS1*NS2 -1);  

  //-------------- ideal gen by C-M matrix is prime?--------------------
  I_M:=ideal<CR2|NewBasis>;
  printf"I_M is Radical: %o\n",IsRadical(I_M);  
end for;

End:=Realtime();
printf"compute time Total=%o\n",End-Start;