/*******************************************************************
*
*   anti-cyclotomic.m   -- test Mazur's "sign conjecture"
*
*   William A. Stein
*
*   April 2001
*
*   The file anti-cyclotomic_height.tex, which Mazur wrote,
*   describes the algorithm that is implemented below.
*
********************************************************************/

declare verbose ac_height, 2;    // since mine isn't in the kernel yet.



intrinsic IsAcceptableK(E::CrvEll, p::RngIntElt, D::RngIntElt) -> BoolElt
//   {True if and only if K=Q(sqrt(D)) satisfies conditions (a)--(d).}
   {True if rank(E) is odd, which for all practical purposes in our cases means 1.}

   if not IsFundamental(D) then
      return false;
   end if;

   // condition (a)
   N := Conductor(E);
   if GCD(D,N*p) ne 1 then     
      return false;
   end if;

   // condition (b)
   R<x> := PolynomialRing(RationalField());
   K := NumberField(x^2-D);
   O := MaximalOrder(K);
   if #Factorization(ideal<O|p>) eq 1 then
      return false;
   end if;


/*
   // condition (c)
   // all primes dividing N split
   for fac in Factorization(N) do
      if IsOdd(fac[2]) and #Factorization(ideal<O|fac[1]>) eq 1 then
         return false;
      end if;
   end for;
*/


   // condition (d)
   E_D := QuadraticTwist(E,D);
print Evaluate(KroneckerCharacter(D),-1);
   if Evaluate(KroneckerCharacter(D),-N) eq 1 then    // there's no chance to have rank 1
      return false;
   end if;

   return true;
end intrinsic;


intrinsic FindAcceptableK(E::CrvEll, p::RngIntElt, Dmin::RngIntElt, Dmax::RngIntElt) -> SeqEnum
   {The fields K whose discriminant lie between Dmin and Dmax and satisfy conditions (a)--(d).}

   require IsPrime(p) and p gt 0 : "Argument 2 must be a prime number.";
   require Dmin lt 0: "Argument 3 must be negative.";
   require Dmax lt 0: "Argument 4 must be negative.";
   require Dmin le Dmax: "Argument 3 must be less than or equal to argument 4.";

   N := Conductor(E);
   vprint ac_height : "Finding acceptable K for E of conductor", N;
   return Reverse([D : D in [Dmin..Dmax] | IsAcceptableK(E,p,D)]);
end intrinsic;


intrinsic FindAcceptableK(E::CrvEll, p::RngIntElt) -> RngIntElt
{}
   D := -1;
   while true do
      if IsAcceptableK(E,p,D) then
         return D;
      end if;
      D := D - 1;
   end while;
end intrinsic;


intrinsic rhotilde_cheat(x::FldNumElt , pi::RngLocElt, pibar::RngLocElt) -> FldLocElt
   {}

   vprint ac_height : "rhotilde_cheat(x,", ",", pi, ",", pibar,")";

   p := Prime(Parent(pi));
   time n := Numerator(Norm(x));
/*   if n mod p eq 0 then 
      vprint ac_height : "Bad point.";
      return false;
   end if;
*/
   O := MaximalOrder(Parent(x));
   u := #UnitGroup(O);
   v := (p-1)*u;
   alpha_h := x^v;
   rep := Eltseq(alpha_h);
   vprint ac_height : "Computing embeddings of canonical power of generator into Qp";
   pi_alpha_v    := (rep[1]+rep[2]*pi)^((p-1)*u);
   pibar_alpha_v := (rep[1]+rep[2]*pibar)^((p-1)*u);

   vprint ac_height : "pi_alpha_v =", pi_alpha_v;
   vprint ac_height : "pibar_alpha_v =", pibar_alpha_v;

   vprint ac_height : "pi_alpha_v-1 = ", pi_alpha_v-1;
   val1 := Valuation(pi_alpha_v-1);
   vprint ac_height : "Valuation(pi_alpha_v-1) = ", val1;
   vprint ac_height : "pibar_alpha_v-1 = ", pibar_alpha_v-1;
   val2 := Valuation(pibar_alpha_v-1);
   vprint ac_height : "Valuation(pibar_alpha_v-1) = ", val2; 
 
   if val1*val2 eq 0 then
      vprint ac_height : "Bad point.";
      return false;
   end if;

   vprint ac_height : "Valuations ok -- now computing p-adic logarithms.";
   logpi := Log(pi_alpha_v);
   vprint ac_height : "lambda(pi_alpha_v) =", logpi;
   logpibar := Log(pibar_alpha_v);
   vprint ac_height : "lambda(pibar_alpha_v) =", logpibar;

   return 1/v * ( logpi - logpibar );
end intrinsic;


intrinsic rhotilde(I::RngOrdIdl , pi::RngLocElt, pibar::RngLocElt) -> FldLocElt
   {I is an ideal of O_K that is coprime to p. Each of pi and pibar define
   a map from O_K to Z_p; pi is the image of O.1, and similarly for pibar.
   (Note that pi and pibar are computed by finding the roots of the charpoly
   of O.1 over Z_p.)}

   vprint ac_height : "rhotilde(I,", pi, ",", pibar,")";

   p := Prime(Parent(pi));
   time n := Numerator(Norm(I));
   if n mod p eq 0 then 
      vprint ac_height : "Bad point.";
      return false;
   end if;
   O := Order(I);
   u := #UnitGroup(O);
   vprint ac_height: "Computing class number h";
   h := ClassNumber(O);
   vprint ac_height: "Class number h =", h;
   v := (p-1)*u*h;
   vprint ac_height : "Raising ideal to the power h";
   Ih  := I^h;
   vprint ac_height : "Finding generator for the principal ideal I^h";

//   alpha_h := Generators(Ih)[2];
   _, alpha_h := IsPrincipal(Ih);
   vprint ac_height : "I^h is principal.";
   rep := Eltseq(alpha_h);
   vprint ac_height : "Computing embeddings of canonical power of generator into Qp";
   pi_alpha_v    := (rep[1]+rep[2]*pi)^((p-1)*u);
   pibar_alpha_v := (rep[1]+rep[2]*pibar)^((p-1)*u);
   vprint ac_height : "pi_alpha_v =", pi_alpha_v;
   vprint ac_height : "pibar_alpha_v =", pibar_alpha_v;

   vprint ac_height : "pi_alpha_v-1 = ", pi_alpha_v-1;
   val1 := Valuation(pi_alpha_v-1);
   vprint ac_height : "Valuation(pi_alpha_v-1) = ", val1;
   vprint ac_height : "pibar_alpha_v-1 = ", pibar_alpha_v-1;
   val2 := Valuation(pibar_alpha_v-1);
   vprint ac_height : "Valuation(pibar_alpha_v-1) = ", val2; 
 
   if val1*val2 eq 0 then
      vprint ac_height : "Bad point.";
      return false;
   end if;

   vprint ac_height : "Valuations ok -- now computing p-adic logarithms.";
   logpi := Log(pi_alpha_v);
   vprint ac_height : "lambda(pi_alpha_v) =", logpi;
   logpibar := Log(pibar_alpha_v);
   vprint ac_height : "lambda(pibar_alpha_v) =", logpibar;

   return 1/v * ( logpi - logpibar );
end intrinsic;

intrinsic Num_and_Denom_Ideals(x::FldNumElt : cheat := false) -> RngOrdIdl, RngOrdIdl
   {I and J, where I is the numerator ideal of x and J is the denominator ideal of x.}
   
   // We compute the denominator ideal as follows.  Compute
   // the intersection of the image of multiplication by x
   // and O inside K.  Let a and b be generators for this intersection,
   // as a Z-module.  Then the denominator ideal is generated by
   // a/x and b/x.  The numerator ideal is then generated by a and b.

   vprint ac_height : "NumeratorIdeal()";

   L   := Parent(x);
   O   := MaximalOrder(L);

   if cheat eq true then  
      d := LCM([Denominator(y) : y in Eltseq(x)]);
      return d*x, L!d;
   end if;

   I   := ideal<O|1>;
   xI  := x*I;
   B   := [L!z : z in Basis(I meet xI)];
   a,b := Explode(B);
   J := a/x*O + b/x*O;
   I := ideal<O|a,b>;
   return I, J;
   
end intrinsic;


intrinsic H_2(P::CrvEllPt, p::RngIntElt : cheat := false) -> FldLocElt
   {The function H_2 defined in Mazur's note.  This is rhotilde(J*x(P))/2, where
   J is the "ideal denominator" of the principal ideal (x(P)) generated by
   the x-coordinate of P.}

   vprint ac_height : "H_2()";
   
   require IsPrime(p) and p gt 0 : "Argument 2 must be a prime number.";
   
   E := Parent(P);
   K := BaseRing(E);
   require Type(K) eq FldNum and Degree(K) eq 2 and Discriminant(K) lt 0 :
         "The base field of the parent of argument 1 must be a quadratic imaginary number field.";

   O := MaximalOrder(K);
   f := PolynomialRing(pAdicRing(p : Precision := 40))!MinimalPolynomial(Basis(O)[2]);
   fac := Factorization(f);
   require #fac eq 2 : "Argument 2 must split in the base field of the parent of argument 1.";
   pi := -Evaluate(fac[1][1],0);
   pibar := -Evaluate(fac[2][1],0);

   if P[3] eq 0 then      
      return pAdicField(p)!0;
   end if;
   x := P[1]/P[3];

   I, J := Num_and_Denom_Ideals(x : cheat := cheat);
   if cheat eq true then
      r := rhotilde_cheat(I,pi,pibar);
   else
      r := rhotilde(I, pi, pibar);
   end if;
   if Type(r) eq BoolElt then
      return r;
   end if;
   return r/2;
   
end intrinsic;


intrinsic Pairing(P::CrvEllPt, Q::CrvEllPt, p::RngIntElt, n::RngIntElt) -> FldLocElt
   {Let h(x) = H_2(p^n*x)/p^(2n).  Then this returns
   h(P+Q) + h(P-Q) - 2*h(P) - 2*h(Q),
   which is an approximation for the height pairing of P and Q.}
   function h(x,p)
      return H_2(p^n*x,p)/p^(2*n);
   end function;
   return h(P+Q,p) + h(P-Q,p) - 2*h(P,p) - 2*h(Q,p);
end intrinsic;


intrinsic PointFromTwist(E::CrvEll, D::RngIntElt) -> CrvEllPt
   {The point on E/Q(sqrt(D)) corresponding to a generator of E_D(Q)/tor,
   assuming that the latter group has rank 1.}

   Ew := WeierstrassModel(E);
   _,_,_,a,b  := Explode(aInvariants(Ew));
   E_D := EllipticCurve([D^2*a,D^3*b]);
   G,f := MordellWeilGroup(E_D : Bound := 2, HeightBound := 0.5);  
   //   require TorsionFreeRank(G) eq 1 : "The twist of argument 1 by argument 2 must have rank 1.";
   if TorsionFreeRank(G) gt 1 then
      print "WARNING: rank is > 1 for E = ", E, " D = ", D;
   end if;
   if TorsionFreeRank(G) lt 1 then
      return false;
      /*
      print "WARNING: rank of E_D computed is < 1.  We have E = ", aInvariants(E), " D = ", D;
      print "WARNING: rank computed is < 1 for E_D = ", aInvariants(E_D);
      G,f := MordellWeilGroup(E_D);  // this might take longer...!
      print "Computed generators and go = ", G;
      */
   end if;

   // P is the generator for E_D(Q)/tor
   P := f(G.Ngens(G));  
   x := P[1]/P[3];
   y := P[2]/P[3];

   R<z> := PolynomialRing(Rationals());
   K<sqrtD> := NumberField(z^2-D);
   Ew,_,psi := WeierstrassModel(E);
   Ew := BaseExtend(Ew,K);
   Pw := [x/D, (y/D^2)*sqrtD];
   Q := Ew!Pw;
   Qx := Q[1]/Q[3];
   Qy := Q[2]/Q[3];   
   r,s,t,u := Explode(IsomorphismData(psi));

   // (x, y) |-> (u^2x + r, u^3 y + su^2 x + t).
   X := u^2*Qx + r;
   Y := u^3*Qy + s*u^2*Qx + t;
   E := BaseExtend(E,K);
   return E![X,Y], E;
   
end intrinsic;

intrinsic TestNontrivialityConjecture(E::CrvEll, p::RngIntElt, 
                                      n::RngIntElt) -> List, RngIntElt
{}
   require IsPrime(p) and p gt 0 : "Argument 2 must be a prime number.";
   require Conductor(E) mod p ne 0 : "Argument 1 must have good reduction at argument 2.";

   D := FindAcceptableK(E,p);
   return TestNontrivialityConjecture(E,p,n,[D]), D;
end intrinsic;


intrinsic TestNontrivialityConjecture(E::CrvEll, p::RngIntElt, n::RngIntElt, Dlist::SeqEnum) 
       -> List
{Let E be a rank two elliptic curve and let P1, P2 be a basis for E(Q).  
 Let P = P1 + P0, where P0 is a point of infinite order in E_D(Q) subset E(Q(sqrt(D))). 
 The value returned is a list of pairs 

          <D, [* H_2(P), H_2(p*P)/p^2, H_2(p^2*P)/p^4, ..., H_2(p^(2*n)*P)/p^(2*n) *]> 

 Some of the H_2 values might equal false, if the point doesn't like in the 
 appropriate subgroup for H_2 to be defined.
}

   require IsPrime(p) and p gt 0 : "Argument 2 must be a prime number.";
   require Conductor(E) mod p ne 0 : "Argument 1 must have good reduction at argument 2.";

   G,f := MordellWeilGroup(E);
   require TorsionFreeRank(G) eq 2 : "Argument 1 must have rank 2.";
   P1 := f(G.(Ngens(G)-1));
   P2 := f(G.Ngens(G));

   ans := [* *];
   for D in Dlist do
      P0 := PointFromTwist(E,D);
      if Type(P0) eq BoolElt then
	 print "Couldn't find a point on the twist quickly enough, so bailing.";
	 print "E = ", E;
	 print "D = ", D;
	 continue;
      end if;
      Pa  := P0 + Parent(P0)!Eltseq(P1);     
      Pb  := P0 + Parent(P0)!Eltseq(P2);     
      hts := [* *];
      for m in [1..n] do
         H2a := H_2(Pa,p);
         H2b := H_2(Pb,p);

         if Type(H2a) ne BoolElt then
            H2a := H2a/p^(2*m);
         end if;
         if Type(H2b) ne BoolElt then
            H2b := H2b/p^(2*m);
         end if;
         Append(~hts, <H2a, H2b>);

         vprint ac_height : "When m =", m, " hts = ", hts;

         if m lt n then
            vprint ac_height : "Multiplying each point by", p;
            time Pa := p*Pa;
            time Pb := p*Pb;
         end if;
      end for;
      Append(~ans, <D, hts>);
      vprint ac_height : "Answer so far = ", ans;
   end for;

   return ans;     
end intrinsic;



intrinsic Rank2Curve(n::RngIntElt) -> CrvEll
{}
   curves := ["389A","433A","446D", "563A","571B","643A","655A","664A","681C",
              "707A","709A","718B","794A","817A","916C","944E","997B","997C","1001C",
              "1028A","1034A","1058C","1070A","1073A","1077A","1088J","1094A","1102A",
              "1126A","1132A","1137A","1141A","1143C","1147A","1171A"];
   require n ge 1 and n le #curves : "Argument 1 must be between 1 and", #curves;

   return EllipticCurve(CremonaDatabase(),curves[n]), curves[n];
end intrinsic;


intrinsic SmallestDTest(p::RngIntElt)
{}
   out := Open(Sprintf("SmallestDTest_p%o.out",p),"w");

   for n in [1..35] do
      E, label := Rank2Curve(n);
      if Conductor(E) mod p eq 0 then
         continue;
      end if;
      f := qEigenform(E,p+1);
      if Coefficient(f,p) eq 0 then
	 continue;
      end if;
      fprintf out, "\n// ** %o **: %o\n", label, E;
      Flush(out);
      ans, D := TestNontrivialityConjecture(E,p,4);
      fprintf out, "\n// Q(sqrt(%o))\n\n", D;
      fprintf out, "heights[%o] := %o;\n\n", n, ans;
      Flush(out);
   end for;
end intrinsic;


intrinsic EquidisTest(curve::MonStgElt, p::RngIntElt, n::RngIntElt, 
             Dstop::RngIntElt, Dstart::RngIntElt : cheat := false)
{}
//   curve := "389A";
   E := EllipticCurve(CremonaDatabase(),curve);
   D := FindAcceptableK(E,p,Dstop,Dstart);
   Qp<q> := pAdicField(p);
   Qp`SeriesPrinting := true;

   G,f := MordellWeilGroup(E);
   out := Open(Sprintf("equidis_test_%o_p%o_n%o_d%o_%o.m", curve, p, n, Dstart,cheat),"w");
   fprintf out, "D = %o\n", D;

   for d in D do
      t := Cputime();
      fprintf out, "\nK = Q(sqrt(%o)):\t", d;
      P  := PointFromTwist(E,d);
      if Type(P) eq BoolElt then 
         continue;
      end if;
      EE := Parent(P);
      P0 := EE!Eltseq(f(G.1));
      P1 := EE!Eltseq(f(G.2));
      hAC := H_2(p^n*(P0+P),p : cheat := cheat);
      if Type(hAC) ne BoolElt then
         hAC := hAC/p^(2*n);
         hBC := H_2(p^n*(P1+P),p : cheat := cheat);
         if Type(hBC) ne BoolElt then
            hBC := hBC/p^(2*n);
         end if;
      end if;
      if Type(hAC) eq BoolElt or Type(hBC) eq BoolElt then
         fprintf out, "(no data)";      
      else
         r := hAC/hBC;
         fprintf out, "%o", r;
      end if;
      fprintf out, " (time = %o)\n", Cputime(t);
      Flush(out);
   end for;
   
end intrinsic;
















/* example usage

Attach("anti-cyclotomic.m");
E := EC("389A"); 
SetVerbose("ac_height",2);
p := 3;

Dp := FindAcceptableK(E,p,-35,-1);
P := [* *];
for d in Dp do
   Append(~P, PointFromTwist(E,d));
end for;
EE := Parent(P[2]);
G,f := MordellWeilGroup(E);
P0 := EE!Eltseq(f(G.1));
P1 := EE!Eltseq(f(G.2));
H_2(P0+P[2],p);
H_2(P1+P[2],p);
*/

