////////////////////////////////////////////////////////////////////////
//                                                                    //
// merel-stein.m -- MAGMA code for verifying conditions 1)-3) in the  //
// paper "The field generated by the points of small prime order on   //
// an elliptic curve."                                                //
//                                                                    //
////////////////////////////////////////////////////////////////////////

// Note: The sequence MOG, LChiZero, and Anomalous are all defined
// in other files available with this file.


intrinsic MerelNumberOfIsogenyClasses(MOG::SeqEnum, p::RngIntElt) -> RngIntElt
{Returns the number of isogeny classes at level p.}
   return #MOG[p]`Eigenvectors;
end intrinsic;

intrinsic MerelNumberOfEmbeddings(MOG::SeqEnum,
                       p::RngIntElt, 
                       IsogenyClass::RngIntElt
                       : Factors := true) -> RngIntElt
{Number of factors mod p of the polynomial defining IsogenyClass i.}
   R := PolynomialRing(GF(p));
   f := R ! MOG[p]`Eigenvectors[IsogenyClass]`DefiningPolynomial;
   if Factors then
      return #Factorization(f);
   end if;
   return &+[Degree(a[1]) : a in Factorization(f)];  
end intrinsic;

intrinsic MerelDimension(MOG::SeqEnum,
                         p::RngIntElt, 
                         IsogenyClass::RngIntElt) -> RngIntElt
{Number of factors mod p of the polynomial defining IsogenyClass i.}
   return Degree(MOG[p]`Eigenvectors[IsogenyClass]`DefiningPolynomial);
end intrinsic;

intrinsic MerelIota(MOG::SeqEnum,
                      p::RngIntElt, 
           IsogenyClass::RngIntElt, 
                      u::RngIntElt,
              Embedding::RngIntElt
                           : Factors := true) -> .
{Given a level p <= 997, an IsogenyClass, an integer u,
 and an Embedding number <= MerelNumberOfEmbeddings(p,IsogenyClass), 
 return the rational function sum n_E/((j-j(E))^u.  This function
 is defined in section 1.4 of Merel's paper "Sur la nature  
 non-cyclotomique des points d'ordre fini des courbes elliptiques."  
 The optional parameter Factors controls how the Embedding number
 is defined.}
   R   := PolynomialRing(GF(p));
   v   := MOG[p]`Eigenvectors[IsogenyClass];
   fac := Factorization(R!v`DefiningPolynomial);
   if Factors then
      f := R!fac[Embedding][1];
      i := 1;
   else
      i   := 1;
      while Embedding gt Degree(fac[i][1]) do
         Embedding -:= Degree(fac[i][1]);
         i         +:= 1;
      end while;
      f   := R!fac[i][1];
      i   := Embedding;
   end if;
   Fq<a>:= GF(p^LCM(2,Degree(f)));    // need the s.s j-invariants.
   RR  := PolynomialRing(Fq);
   alp := Roots(RR!f)[i][1];            // assuming output of Roots sorted.
   ss_f:= MOG[p]`SupersingularBasis`DefiningPolynomial;
   bet := Roots(RR!ss_f)[1][1];         // first is good enough.
   coor:=MOG[p]`SupersingularBasis`Coordinates;
   ss_j:= [ Evaluate(RR!coor[n],bet) : n in [1..#coor] ];
   FF<j>:= FieldOfFractions(RR);   
   return &+ [
      FF!Evaluate(RR!v`Coordinates[n],alp)           // n_e
       / (FF!ss_j[n]-j)^u  : n in [1..#ss_j]
   ];
end intrinsic;

////////////////////////////////////////////////////////////
// Criterion when p = 1 (mod 4):
// It turns out that this more complicated criterion, which
// was mentioned in Merel's first paper [2], is not needed.
// We instead only must verify the criterion 
// "for p=3 (mod 4)" below.
////////////////////////////////////////////////////////////
function MerelCriterion1Degree(MOG, LChiZero, Anomalous, p, d)
   assert p mod 4 eq 1;
   assert d gt 0 and (p-1) mod d eq 0;
   assert d ne 2;

   NumberOfClasses := MerelNumberOfIsogenyClasses(MOG,p);
   Anom            := Anomalous[p];
   for iso1 in [1..NumberOfClasses] do
      for iso2 in [iso1..NumberOfClasses] do
         if Index(LChiZero[p][iso1],d) eq 0 and
            Index(LChiZero[p][iso2],d) eq 0 then  
            for embedding1 in 
              [1..MerelNumberOfEmbeddings(MOG,p,iso1 : Factors:=false)] do
               for embedding2 in 
                  [1..MerelNumberOfEmbeddings(MOG,p,iso2 : Factors:=false)] do
                  f11   := MerelIota(MOG, p, iso1, 1, embedding1 : Factors := false);
                  f12   := MerelIota(MOG, p, iso1, 2, embedding1  : Factors := false);
                  f21   := MerelIota(MOG, p, iso2, 1, embedding2 : Factors := false);
                  f22   := MerelIota(MOG, p, iso2, 2, embedding2 : Factors := false);
                  // Magma is clever enough to force the following to
                  // to take place in a common finite extension of
                  // the finite fields over which each fij is defined. 
                  f     := f11*f22 - f12*f21; 
//"f11=",f11;
//"f12=",f12;
//"f21=",f21;
//"f22=",f22;
//"f=",f;
                  vals  := [Evaluate(f,j) : j in Anom];
                  if Index(vals,0) eq 0 then
                     printf "verified with A=%o%o(dim=%o), B=%o%o(dim=%o)\n",  
                       p, ToIsogenyCode(iso1), MerelDimension(MOG,p,iso1),
                       p, ToIsogenyCode(iso2), MerelDimension(MOG,p,iso2);
                     return true;
                  end if;
               end for;
            end for;
         end if;
      end for;
   end for;
   return false;
end function;

intrinsic MerelCriterion1(MOG::SeqEnum, LChiZero::SeqEnum,
                          Anomalous::SeqEnum,  p::RngIntElt
                          ) -> BoolElt
{}
   assert p mod 4 eq 1;
   "--------->>>>>>  p =",p, " <<<<<<<---------";
   for d in [e : e in Divisors(p-1) | e ne 2] do
      printf "deg(chi)=%3o   ",d;
      if not MerelCriterion1Degree(MOG, LChiZero, Anomalous,p,d) then
         "FAILED!";
         return false;
      end if;
   end for;
   "TRUE.";
   return true;
end intrinsic;


////////////////////////////////////////////////////////////
// Criterion when p = 3 (mod 4):
function MerelCriterion3Degree(MOG, LChiZero, Anomalous, p, d)
   /* assert p mod 4 eq 3; */   // commented out, because this
                                // criterion is also valid for p=1(mod 4).
   assert d gt 0 and (p-1) mod d eq 0;
   NumberOfClasses := MerelNumberOfIsogenyClasses(MOG,p);
   Anom            := Anomalous[p];
   for iso in [1..NumberOfClasses] do
      if Index(LChiZero[p][iso],d) eq 0 then   // L(iso,[d],1) =/= 0.
         for embedding in [1..MerelNumberOfEmbeddings(MOG,p,iso)] do
            f    := MerelIota(MOG, p, iso, 1, embedding);
            vals := [Evaluate(f,j) : j in Anom];
            if Index(vals,0) eq 0 then
               printf "verified with %o%o, (dim=%o)\n",       
                       p, ToIsogenyCode(iso), 
                       MerelDimension(MOG,p,iso);
               return true;
            end if;
         end for;
      end if;
   end for;
   return false;
end function;

intrinsic MerelCriterion3(MOG::SeqEnum, LChiZero::SeqEnum,
                               Anomalous::SeqEnum,  p::RngIntElt
                         ) -> BoolElt
{}
   /* assert p mod 4 eq 3; */
   "--------->>>>>>  p =",p, " <<<<<<<---------";
   for d in Divisors(p-1) do
      printf "deg(chi)=%3o   ",d;
      if p mod 4 eq 1 and d eq 2 then
         "Quadratic, hence skipping...";
      else
         if not MerelCriterion3Degree(MOG, LChiZero, Anomalous,p,d) then
            "FAILED!";
            return false;
         end if;
      end if;
   end for;
   "TRUE.";
   return true;
end intrinsic;

intrinsic MerelCriterion(MOG::SeqEnum, LChiZero::SeqEnum,
                               Anomalous::SeqEnum,  p::RngIntElt) -> BoolElt
{Verifies the criterion of Merel for the prime number p.
 The sequence MOG, LChiZero, and Anomalous are all defined
 in other files available with this file.}
   assert IsPrime(p) and p ge 11;
   return MerelCriterion3(MOG,LChiZero,Anomalous,p);
   /*  because the 3 mod 4 criterion suffices for both cases, as it turns out. 
   case p mod 4:
      when 1: return MerelCriterion1(MOG,LChiZero,Anomalous,p);
      when 3: return MerelCriterion3(MOG,LChiZero,Anomalous,p);
   end case;
   */
end intrinsic;