/* bsd.m -- MAGMA Implementation of algorithm to compute Sha(E/Q)[p]
   for all primes p for which there are no p-isogenies. */

function IntSeq(P)
   assert Type(P) eq PtEll;
   if P[3] eq 0 then
      return [0,1,0];
   end if;
   d := LCM([Denominator(P[i]) : i in [1..3]]);
   return [Integers() | P[1]*d, P[2]*d, P[3]*d];
end function;

intrinsic ShaPrimes(E::CrvEll, D::RngIntElt) -> RngIntElt
{Returns the product of primes that might divide Sha(E/Q).  Returns
0 to indicate no information.}
   print "Using Kolyvagin theorem on D =", D;
   tm := Cputime();
   E := MinimalModel(E);
   print "aInvariantsE=", aInvariants(E);
   K := QuadraticField(D);
   OK := MaximalOrder(K);
   u := 1;   // Half the number of units of OK.
   D := Discriminant(K);
   require not (D in [-4, -3]):  "Extra units not allowed.";
  
   // this is a technical hypothesis needed (by us, at least!) to
   // imply that hypothesis to Kolyvagin are satisfied. 
   //require #Factorization(D) ge 2 : "Discriminant must be divisible by 2 primes.";

   rE, LE := AnalyticRank(E);
   require rE le 1: "Analytic rank of E must be at most 1.";
   print "rE=", rE;
   print "LE=", LE;

   F := MinimalModel(QuadraticTwist(E,D));
   print "aInvariantsF=", aInvariants(F);
   rF, LF := AnalyticRank(F);
   print "rF=", rF;
   print "LF=", LF;
   require rF le 1: "Analytic rank of twist must be at most 1.";

   require rE + rF eq 1: "Sum of analytic ranks must be 1.";

   N := Conductor(E);
   print "N = ", N;
   require GCD(D, N) eq 1: "D must be coprime to conductor of E.";

   for p in Factorization(N) do
      require #Factorization(p[1]*OK) eq 2 :
            "Every prime dividing conductor must split.";
   end for;


   B := 1;

   // 0. If D is divisible by only one prime, exclude it.
   if #Factorization(D) eq 1 then
      print "D =",D," is divisible by only prime, so excluding that prime from B..";
      B *:= Factorization(D)[1][1];
   end if;

   // 1. Find all curves in the isogeny class and exclude order
   // of torsion subgroup.

   I := IsogenousCurves(E);
   printf "IsogenousCurves=[";
   for i in [1..#I] do
      printf "%o",aInvariants(I[i]);
      if i ne #I then
         printf ",";
      end if;
   end for;
   printf "]\n";
   B := LCM(2,B);  // always exclude 2.
   torK := [#TorsionSubgroup(BaseExtend(F,K)) : F in I];
   print "torK=",torK;
   for t in torK do
      B := LCM(B, t);
   end for;

   // 3. Root number
   eps := RootNumber(E);   // "epsilon"
   print "eps_E=", eps;

   /* 4. [Compute Mordell-Weil] If the root number is $-1$, compute
      $E(\Q)$, and let $z$ be a generator modulo torsion.  Otherwise, 
      if the root number is $+1$, compute $E^D(\Q)$, and let $z$ be a
      generator modulo torsion. */
  if eps eq -1 then
     G, f := MordellWeilGroup(E);
     if Rank(E) ne 1 then
        print "WARNING: MAGMA definitely failed to compute rank of E correctly!!!!";
        return 0, aInvariants(F), [0,0,0], 0;
     end if;
  else
     G, f := MordellWeilGroup(F);
     if Rank(F) ne 1 then
        print "WARNING: MAGMA definitely failed to compute rank of F correctly!!!!";
        return 0, aInvariants(F), [0,0,0], 0;
     end if;
  end if;
  z := f(G.Ngens(G));
  print "MordellWeilF=",IntSeq(z);

  P := Periods(E);
  omega := 2*Abs( Real(P[1])*Imaginary(P[2]) - Imaginary(P[1])*Real(P[2]) );
  print "omega_E=",omega;
  
  // 5. Height of the Heegner Point.
  alpha := u^2 * Sqrt(Abs(D)) / omega;
  height_heegner_over_K := alpha * LE * LF;
  
  // 6. Index of Heegner point.
  height_of_z_over_K := 2 * Height(z);
  a := height_heegner_over_K / height_of_z_over_K;

  a := a*4;
  
  b := Round(a);
  print "Rounding... error in index computation: ", Abs(b-a);

  // Note -- we have to worry about precisions of computations, rounding
  // errors, etc., to make this all rigorous.

  it_is_a_square, sqrt_b := IsSquare(b);
  print "quotient = ", a;
  if not it_is_a_square then
     print "Because of rounding errors, trying +1 and -1 in extracting square root when computing index.";
     // Try b-1 and b+1:
     it_is_a_square, sqrt_b := IsSquare(b+1); 
     if not it_is_a_square then
        it_is_a_square, sqrt_b := IsSquare(b-1); 
     end if;
     if not it_is_a_square then
        print "The quotient h_heeg / Height(z) is ", a, " which rounds to", b,", which is not a square, or within 1 of a square!"; 
        print "WARNING!!!";
        return 0, aInvariants(F), [0,0,0], 0;
     end if;
  end if;
  I_K := sqrt_b;
  print "Index=", I_K;
  print "B= ", B;
  if I_K eq 0 then
     print "WARNING -- index of 0!!";
     print "The quotient h_heeg / Height(z) is ", a, " which rounds to a number with sqrt 0."     ;
     return 0, aInvariants(F), [0,0,0], 0;
  end if;

  bound := &*[p[1] : p in Factorization(I_K*B)];
  print "Obtained the bound ", bound, " in ", Cputime(tm), " seconds.";
  return bound, aInvariants(F), IntSeq(z), Height(z);

end intrinsic;


intrinsic ShaPrimes(E::CrvEll, D_range::SeqEnum) -> SeqEnum
{Return bound on primes that divide Sha using all
D's in D_range.}

   N := Conductor(E);
   good_disc := [];
   for D in D_range do
      if not IsFundamentalDiscriminant(D) then
         continue;
      end if;
      K := QuadraticField(D);
      OK := MaximalOrder(K);
      disc := Discriminant(K);
      if #Factorization(disc) le 1 then
         print disc, " fails because it is not divisible by 2 primes.";
         continue;
      end if;
      if disc in [-4, -3] then
         print disc, " fails Heegner since is -3, -4.";
         continue;
      end if;
      if GCD(disc, N) ne 1 then
         print disc, " fails Heegner since is not coprime to conductor.";
         continue;
      end if;
      F := MinimalModel(QuadraticTwist(E,disc));
      rF, LF := AnalyticRank(F);
      if rF ge 2 then
         print disc, " fails Heegner since rank of twist is too big.";
         continue;
      end if;
      splits := true;
      for p in Factorization(N) do
         if #Factorization(p[1]*OK) ne 2 then
            print disc, " fails Heegner since not every prime dividing conductor splits";
            splits := false;
            continue;
         end if;
      end for;
      if splits then
         Append(~good_disc, disc);
      end if;
   end for;

   if #good_disc eq 0 then
      return 0, [];
   end if;
   Bounds := [<ShaPrimes(E,D),D> : D in good_disc];
   return GCD([b[1] : b in Bounds]), Bounds;
end intrinsic;

intrinsic ShaPrimes(E::CrvEll) -> RngIntElt, RngIntElt
{Return bound on primes that divide Sha using first Heegner D.
Also returns the first Heegner D divisible by two primes.}

   N := Conductor(E);
   D := -1;
   while true do
      D -:= 1;
      if not IsFundamentalDiscriminant(D) then
         continue;
      end if;
      K := QuadraticField(D);
      OK := MaximalOrder(K);
      disc := Discriminant(K);
      if disc in [-4, -3] then
         print disc, " fails Heegner since is -3, -4.";
         continue;
      end if;
      if #Factorization(disc) le 1 then
         print disc, " fails because only divisible by 1 prime.";
         continue;
      end if;
      if GCD(disc, N) ne 1 then
         print disc, " fails Heegner since is not coprime to conductor.";
         continue;
      end if;
      F := MinimalModel(QuadraticTwist(E,disc));
      rF, LF := AnalyticRank(F);
      if rF ge 2 then
         print disc, " fails Heegner since rank of twist is too big.";
         continue;
      end if;
      splits := true;
      for p in Factorization(N) do
         if #Factorization(p[1]*OK) ne 2 then
            print disc, " fails Heegner since not every prime dividing conductor splits";
            splits := false;
            continue;
         end if;
      end for;
      if splits then
         return D, ShaPrimes(E, D);
      end if;
   end while;
end intrinsic;

intrinsic ShaPrimes2(E::CrvEll : num_allowed_failures := 1000) -> RngIntElt, RngIntElt, 
RngIntElt, RngIntElt, RngIntElt

{Return bounds on primes that divide Sha using first two Heegner
discriminants D1 and D2.  Also returns D1 and D2.  Thus a bound is the
GCD of the first two returned arguments.

  Return values:
   B -- gcd(B1,B2)
   B1 -- bound from discriminant D1
   B2 -- bound from discriminant D2
   D1 -- discriminant D1
   D2 -- discriminant D2
   z1 -- presumed generator for full MW group for D1
   h1 -- Neron-Tate canonical height of z1
   z2 -- presumed generator for full MW group for D2
   h2 -- Neron-Tate canonical height of z2

  If D1 is divisible by at least two primes, then D2=D1, B2=B1, and
  B=B1.
}

   N := Conductor(E);
   D := -1;
   B1 := -1;
   num_failures := 0;
   while true do
      D -:= 1;
      if not IsFundamentalDiscriminant(D) then
         continue;
      end if;
      K := QuadraticField(D);
      OK := MaximalOrder(K);
      disc := Discriminant(K);
      assert disc eq D;
      if disc in [-4, -3] then
         print disc, " fails Heegner since is -3 or -4.";
         continue;
      end if;
      if GCD(disc, N) ne 1 then
         print disc, " fails Heegner since is not coprime to conductor.";
         continue;
      end if;
      F := MinimalModel(QuadraticTwist(E,disc));
      rF, LF := AnalyticRank(F);
      if rF ge 2 then
         print disc, " fails Heegner since rank of twist is too big.";
         continue;
      end if;
      splits := true;
      for p in Factorization(N) do
         if #Factorization(p[1]*OK) ne 2 then
            print disc, " fails Heegner since not every prime dividing conductor splits";
            splits := false;
            continue;
         end if;
      end for;
      if splits then
         if B1 eq -1 then
            B1,F1,z1,h1 := ShaPrimes(E, D);       
            if B1 eq 0 and num_failures lt num_allowed_failures then
               num_failures +:= 1;
               print "Failed using D = ", D;
               B1 := -1;
               continue;
            end if;
            D1 := D;
            if #Factorization(D) gt 1 then
               B2 := B1;
               D2 := D1;
               z2 := z1;
               h2 := h1;
               F2 := F1;
               break;
            end if;
         else 
            B2,F2,z2,h2 := ShaPrimes(E, D);       
            if B2 eq 0 and num_failures lt num_allowed_failures then
               num_failures +:= 1;
               print "Failed using D = ", D;
               B2 := -1;
               continue;
            end if;
            D2 := D;
            break;
         end if;
      end if;
   end while;
   return GCD(B1,B2), B1, B2, D1, D2, F1, z1, h1, F2, z2, h2;
end intrinsic;

intrinsic AllCurves(Nstart::RngIntElt, Nstop::RngIntElt) 
{}
   db := CremonaDatabase();
   file := Open(Sprintf("shaprimes-%o-%o.data",Nstart,Nstop),"a");
   for N in [Nstart..Nstop] do
      for i in [1..NumberOfIsogenyClasses(db,N)] do
         E := EllipticCurve(db, N, i, 1);
         bound, disc := ShaPrimes(E);
         print "**************************************";
         print "** N =", N, "      i =", i;
         print "** bound =", bound, "    disc =", disc;
         print "**************************************";       
         fprintf file, "%o,\n", [N, i, bound, disc];
      end for;
   end for;
end intrinsic;