/**********************************************************************************
  heegner.m  -- Heegner points package for MAGMA.

  William Stein, June 2002.

  Version 0.1

 Remark: 
  12/18/02 -- fixed some things so it runs in MAGMA V2.9-25

 **********************************************************************************/




/* NOTES:
   * Reminder -- to increase real number precision:  AssertAttribute(FldPr,"Precision",50);

   * Examples:

> E := 

> E := EllipticCurve([0,0,1,-1,0]);   // 37A
> A := [-D : D in [1..100] | HeegnerHypothesis(E,-D)]; A;
[ -3, -4, -7, -11, -40, -47, -67, -71, -83, -84, -95 ]
P := [* *]; for D in A do time Y := HeegnerPoint(E,D); Append(~P,Y); end for;
Time: 0.040
Time: 0.030
Time: 0.040
Time: 0.050
Time: 0.090
Time: 0.180
Time: 0.050
Time: 0.250
Time: 0.120
Time: 0.140
WARNING: Precision overflow or point is 0.  Use L-series to decide or try 
AssertAttribute(FldPr,"Precision",n) or tweaking parameters.
Time: 0.270

> for i in [1..10] do Recognize(P[i]); end for;
(-1 : 0 : 1)
(1 : -1 : 1)
(0 : 0 : 1)
(0 : -1 : 1)
(1 : -1 : 1)
(0 : 0 : 1)
(6 : -15 : 1)
(0 : -1 : 1)
(0 : 0 : 1)
(0 : 0 : 1)


* BUG!?

> E := EllipticCurve([ 0, 1, 1, -26, -61 ]);
> time P := HeegnerPoint(E,-107); P;
> (3.665647323987715053919669397 + 6.116738926164336860764880743 E-25*i : -0.49999999999999999999999867820523179988 + 9.6625393082244036119397240285848446462*i : 1)
> IdentifyPoint(P);
> false 6938055/1892723
// This time it doesn't identify it as a rational point.  In fact, 
// the coordinates of P generate the Hilbert class field of Q(sqrt(-107)),
// which is really weird because we just took a trace to get down
// from the Hilbert class field.  What is going on!  
// Note: Pete Green's Heegner point package gives the same wrong answer.

*/

declare verbose Heegner, 2;
declare attributes CrvEll:
         qexp, real_qexp_over_n,
         wp, wpp, wp_curve, wp_map, wp_prec;

intrinsic qExpansion(E::CrvEll, prec::RngIntElt) -> RngSerPowElt
{q-expansion}
   R<q> := PowerSeriesRing(Integers());
   if not assigned E`qexp then
      E`qexp := q + O(q^2);
      E`real_qexp_over_n := [Exp(0)];
   end if;
   if prec gt AbsolutePrecision(E`qexp) then
      vprintf Heegner: "Computing q-expansion (prec =%o)\n", prec;
      F := MinimalModel(E);
      D := Integers()!Discriminant(F);
      f := &+[Coefficient(E`qexp,n)*q^n : n in [1..AbsolutePrecision(E`qexp)-1]];

      aplist := TracesOfFrobenius(E, Max(2,prec-1));
      prime_pos := 1;
      p := 2;
      while p lt AbsolutePrecision(E`qexp) do
         prime_pos := prime_pos+1;
         p := NextPrime(p);
      end while;

      for n in [AbsolutePrecision(E`qexp)..prec-1] do
         if n eq 1 then
            an := 1;
         elif IsPrime(n) then
            an := aplist[prime_pos];
            prime_pos := prime_pos + 1;
         else
            fac := Factorization(n); 
            if #fac eq 1 then
               p  := fac[1][1];
               r  := fac[1][2];
               an := Coefficient(f,p) * Coefficient(f,p^(r-1))
                        - (D mod p eq 0 select 0 else p*Coefficient(f,p^(r-2)));
            else  // a_m*a_r := a_{mr} and we know all a_i for i<n.
               m  := fac[1][1]^fac[1][2];
               an := Coefficient(f,m)*Coefficient(f,n div m);
            end if;            
         end if;
         f +:= an * q^n;
      end for;
      E`qexp := f + O(q^prec);
      E`real_qexp_over_n := [RealField()!Coefficient(f,n)/n : n in [1..prec-1]];
      vprint Heegner :  "done computing q-expansion";
   end if;
   return E`qexp + O(q^prec), [E`real_qexp_over_n[i] : i in [1..prec-1]];
end intrinsic;

intrinsic LOne(E::CrvEll : prec := -1) -> FldPrElt, RngIntElt
{The value L(E,1) computed using prec terms of q-expansion and the sign of the functional equation.}
   C := ComplexField();
   pi := Pi(C);
   N := Conductor(E);
   if prec eq -1 then
      prec := N+10;
   end if;
   e := Exp(-2*pi/Sqrt(N));
   _, an_over_n := qExpansion(E, prec);
   V := VectorSpace(ComplexField(),prec-1);
   if N eq SquareFree(N) then
      eps := -&*[-Coefficient(qExpansion(E,P[1]+1),P[1]) : P in Factorization(N)];
      if eps eq -1 then
         return C!0, eps;
      end if;
   else
      e1 := Exp(-2*pi*(1.1)/Sqrt(N));
      e2 := Exp(-2*pi/((1.1)*Sqrt(N)));
      L1 := InnerProduct(V!an_over_n, V![e1^n : n in [1..prec-1]]);
      L2 := InnerProduct(V!an_over_n, V![e2^n : n in [1..prec-1]]);
      if Abs(L1-L2) lt 10^(-4) then
         eps := -1;
         return C!0, eps;
      end if;
   end if;
   ans := 2*InnerProduct(V!an_over_n, V![e^n : n in [1..prec-1]]);
   eps := +1;
   if Abs(ans) lt 10^(-10) then
      ans := C!0;
   end if;
   return ans, eps;
end intrinsic;

intrinsic LPrimeOne(E::CrvEll : prec := -1) -> FldPrElt
{The value L'(E,1).}
   N := Conductor(E);
   if prec eq -1 then
      prec := N+10;
   end if;
   L, eps := LOne(E : prec := prec);
   C := ComplexField();
   if eps eq +1 then
      return C!0;
   end if;
   _, an_over_n := qExpansion(E, prec);
   V := VectorSpace(ComplexField(),prec-1);
   alpha := 2*Pi(C)/Sqrt(N);
   ans := 2*InnerProduct(V!an_over_n, V![ExponentialIntegralE1(alpha*n) : n in [1..prec-1]]);
   if Abs(ans) lt 10^(-10) then
      ans := C!0;
   end if;
   return ans;
end intrinsic;


intrinsic LOneChi(E::CrvEll, D::RngIntElt : prec := -1) -> FldPrElt
{}
   return LOne(MinimalModel(QuadraticTwist(E,D)) : prec := prec);
/* 
   N := Conductor(E);
   require GCD(D, N) eq 1 : "Argument 2 and the conductor of argument 1 must be coprime.";
   C := ComplexField();
   L, eps := LOneChi(E);
   if KroneckerSymbol(D,-N) eq eps then
      return C!0;
   end if;
*/
/*
   e := Exp(-2*pi/(D*Sqrt(N)));
   _, an_over_n := qExpansion(E, prec);
   V := VectorSpace(ComplexField(),prec-1);
*/
   
end intrinsic;

intrinsic LPrimeOneChi(E::CrvEll, D::RngIntElt : prec := -1) -> FldPrElt
{}
   return LPrimeOne(MinimalModel(QuadraticTwist(E,D)) : prec := prec);
end intrinsic;


intrinsic WeierstrassP(E::CrvEll, prec::RngIntElt) -> FldPrElt
{}
   vprint Heegner :  "Computing WeierstrassP";
   L<w> := PolynomialRing(Rationals());
   R<t> := LaurentSeriesRing(L);
   p := t^(-2);
   c4, c6 := Explode(cInvariants(E));
   a := 4; b := -c4/12; c := -c6/216;
   for n in [1..prec] do
      ptest := p + w*t^n + O(t^(n+1));
      ptest_prime := Derivative(ptest);
      f := ptest_prime^2 - a*ptest^3 - b*ptest - c;   
      x := Coefficient(f,-4+n);
      coeff := -Coefficient(x,0)/Coefficient(x,1);
      p := p + coeff*t^n;
   end for;
   vprint Heegner :  "Done computing Weierstrass P";
   return LaurentSeriesRing(Rationals())!(p + O(t^(prec+1)));
end intrinsic;

intrinsic z_to_point(E::CrvEll, z::FldPrElt, wp_prec::RngIntElt  : t := 6) -> SeqEnum
{}
   require E eq MinimalModel(E) : "Argument 1 must be minimal.";
   if not assigned E`wp or wp_prec gt E`wp_prec then   
      E`wp := WeierstrassP(E, wp_prec);
      E`wpp := Derivative(E`wp);
      E`wp_prec := wp_prec;

      c4,c6 := Explode(cInvariants(E));
      F := EllipticCurve([-c4/3,-2*c6/27]);
      _,phi,_ := MinimalModel(F);
      G := EllipticCurve([ComplexField()|a : a in aInvariants(F)]);
      E`wp_curve := G;
      E`wp_map := IsomorphismData(phi);
   end if;
   point := [4*Evaluate(E`wp,z/2^t), 4*Evaluate(E`wpp,z/2^t)];
   vprint Heegner: "Point = ", point;
   G := E`wp_curve;
   pointC := 2^t*(G!point);
   // (x,y) |---> (u^2x + r, u^3y + su^2x + t).
   r,s,t,u := Explode(E`wp_map);
   x,y := Explode(Eltseq(pointC));
   imagepoint := [u^2*x + r, u^3*y + s*u^2*x + t];
   return imagepoint;
end intrinsic;

function CoerceIntoC(x)
   _,prec := HasAttribute(FldPr, "Precision");
   return ComplexField()!(ComplexField(prec)!x);
end function;

intrinsic tau_to_z(E::CrvEll, tau::FldQuadElt, qexp_prec::RngIntElt) -> FldPrElt
{}
   C<i> := ComplexField();
   tau := CoerceIntoC(tau);
   if Imaginary(tau) lt 0 then
      tau := ComplexConjugate(tau);
   end if;
   Pi := Pi(C);
   f,an_over_n := qExpansion(E,qexp_prec);
   e := Exp(2*Pi*i*tau);
   V := VectorSpace(C,qexp_prec-1);
   z := InnerProduct(V!an_over_n, V![e^n : n in [1..qexp_prec-1]]);
   return z;
end intrinsic;

intrinsic ideal_to_tau(a::RngQuadIdl, NN::RngQuadIdl, K::FldQuad) -> FldQuadElt
{Given a nonzero integral fractional ideal a of O_K and an integral 
ideal NN such that O_K/NN = Z/NZ, compute tau in K such that
(E(tau),C(tau)) = (C/a^(-1), N^(-1)*a^(-1)/a^(-1)).}
   // reduce ideal: gives equivalent point, since just 
   // multiplies a be a scalar.
    a := Ideal(Reduction(QuadraticForm(a)));

   // compute I1 and I2 in point (I1, I2/I1)
   I1 := a^(-1);
   I2 := I1*NN^(-1);

   // Find a basis for I2 of the form omega_1, omega_2/N for
   // some basis om1 and om2 of I1.   Such a basis
   // exists because I2/I1 = Z/NZ.  We compute it using 
   // Smith Normal Form.
   V := VectorSpace(RationalField(),2);
   W := VectorSpaceWithBasis([V | Eltseq(b) : b in Basis(I2)]);
   B1 := Basis(I1);
   A := MatrixAlgebra(Integers(),2)!&cat[Coordinates(W,V!Eltseq(b)) : b in B1];

   // The smith form S of A, together with unimodular 
   // matrices P and Q such that P * A * Q = S
   S, P, Q := SmithForm(A);
   // Apply P to the basis B1 to get om1, om2.
   om1 := P[1,1]*B1[1] + P[1,2]*B1[2];
   om2 := P[2,1]*B1[1] + P[2,2]*B1[2]; 

   // Finally, the point z is 1/om2;
   z := K!(1/om2);
   return z;
end intrinsic;

intrinsic ideal_to_tau(f::QuadBinElt) -> FldQuadElt
{}
    D := Discriminant(f);
    K := QuadraticField(D);
    sqD := K.1;
    if D mod 4 eq 0 then
       sqD := 2*sqD;
    end if;
    z := (-f[2] + sqD)/(2*f[1]);
    return z;
end intrinsic;

intrinsic levelNforms(N::RngIntElt, D::RngIntElt) -> SeqEnum
{As in Step 11 of Stephens's paper.}
    beta := Integers()!Sqrt(Integers(4*N)!D);
    Q := QuadraticForms(D);
    G,m := ClassGroup(Q);
    F := Sort([Reduction(m(g)) : g in G]);
    levelforms := [Q!1 : i in [1..#F]];
    known     := [false : i in [1..#F]];
    found := 0;
    A := 1;
    N4 := 4*N;
    alpha := 0;
    while found lt #F do
       B := beta + N4*alpha;
       C := (B^2 - D)/(N4*A);
       if Denominator(C) eq 1 then
          test := Q![Integers()|A*N,B,C];
          i := Index(F, Reduction(test));
          if not known[i] then
             found := found + 1;
             known[i] := true;
             levelforms[i] := test;
          end if;
       end if;
       alpha := alpha + 1;
       if alpha gt A+2 then
          alpha := 0;
          A := A + 1;
       end if;
    end while;

    return levelforms;
end intrinsic;

intrinsic HeegnerPointUseIdeals(E::CrvEll, D::RngIntElt,
                  q_prec::RngIntElt, wp_prec::RngIntElt) -> SeqEnum, CrvEll
{}
   require IsFundamentalDiscriminant(D) : "Argument 2 must be a fundamental discriminant.";
   require D lt 0 : "Argument 2 must be negative.";
   N := Conductor(E);
   require GCD(D, 4*N) eq 1 : "D and 4*N must be coprime";
   require Type(BaseRing(E)) eq FldRat : "Argument 1 must be defined over the rational field.";

   K := QuadraticField(D); 
   OK := MaximalOrder(K);
   // Find an ideal NN with O/NN = Z/NZ.
   NN := 1*OK;
   for p in Factorization(N) do 
      Q := Factorization(p[1]*OK)[1][1];
      if Degree(Q) ne 1 then
         require false : "Each prime dividing argument 3 must split completely in Q(sqrt(D)).";
      end if;
      NN := NN * Q^p[2];
   end for;
   G, m := ClassGroup(OK);
   vprint Heegner :  "Class number =", #G;
   taus := [ideal_to_tau(m(g), NN, K) : g in G];
   vprint Heegner :  "taus =", taus;
   // taus2 := [ideal_to_tau_2(m(g), NN, K) : g in G];
   // vprint Heegner :  "taus2 =", taus2;
   z := &+[tau_to_z(E, tau, q_prec) : tau in taus];
   vprint Heegner :  "z    =", z;
   // z2 := &+[tau_to_z(E, tau, q_prec) : tau in taus2];
   // vprint Heegner :  "z2    =", z2;
   p := z_to_point(E, z, wp_prec) ;
   vprint Heegner :  "p = ", p;
   // p2 := z_to_point(E, z2, wp_prec) ;
   // vprint Heegner :  "p2 = ", p2;
   return p;
end intrinsic;

intrinsic HeegnerHypothesis(E::CrvEll, D::RngIntElt) -> BoolElt
{True if D satisfies the Heegner Hypotheses, so there is a Heegner 
point attached to discriminant D.} 
   if not IsFundamentalDiscriminant(D) or D ge 0 then
      return false;
   end if;
   N := Conductor(E);   
   return IsSquare(Integers(4*N)!D) and GCD(N,D) eq 1;
end intrinsic;

intrinsic RootNumber(E::CrvEll) -> RngIntElt
{The sign in the functional equation for E, so the analytic rank of E is even
if and only if the root number is +1.}
   require Type(BaseRing(E)) eq FldRat : "Argument 1 must be defined over the rational field.";
   eps := -1;
/*
   for p in Factorization(Conductor(E)) do 
   end for;
*/
end intrinsic;

intrinsic HeegnerPoint(E::CrvEll, D::RngIntElt :
                        q_prec:= -1, wp_prec := 40, t := 6) -> PtEll, PtEll
{The Heegner point on E corresponding to D (over the complex numbers), and the a corresponding Q-rational point on 
 either E or on the twist E^D (this point might be a trace from Q(sqrt(D)) to Q).  Returns 0 if the Heegner point is torsion.}
   require E eq MinimalModel(E) : "Argument 1 must be minimal.";
   require HeegnerHypothesis(E,D) : "The Heegner Hypotheses must hold.";
   EE := EllipticCurve([ComplexField()|a : a in aInvariants(E)]);
   tiny := 10^(-6);

   /*
   if Abs(LOne(E)) lt tiny then  // E has rank > 0
      if Abs(LPrimeOne(E)) lt tiny or Abs(LOneChi(E,D)) lt tiny then
         return EE!0, E!0; 
      end if;
   else      // E has rank 0, so need twist to have rank at most 1
      if LPrimeOneChi(E,D) lt tiny then
         return EE!0, E!0; 
      end if;
   end if;   
   */

   N := Conductor(E);
   require Type(BaseRing(E)) eq FldRat : "Argument 1 must be defined over the rational field.";

   if q_prec le 0 then
      q_prec := 2*N*Ceiling(Log(-D)/Log(5)) + 100;
      vprint Heegner : "Using q-expansion to precision", q_prec;
   end if;

   forms := levelNforms(N, D);
   taus := [ideal_to_tau(f) : f in forms];
   vprint Heegner :  "h_K   =", #taus;
   vprint Heegner :  "tau's =", taus;
   z := &+[tau_to_z(E, tau, q_prec) : tau in taus];
   vprint Heegner :  "z     =", z;
   p := z_to_point(E, z, wp_prec : t := t) ;
   vprint Heegner :  "p     =", p;
   if not IsPoint(EE,p) then
      error "*WARNING*: Precision Overflow!";
   end if;
   p := EE!p;

   if Abs(LOne(E)) lt tiny then  // E has rank >=1
      if Max(Abs(Imaginary(p[1])), Abs(Imaginary(p[2]))) gt 10^(-2) then
         q := TracePoint(p);
      else
         q := p;
      end if;
      t, Q := Recognize(q);
   else   // E has rank 0 
      phi := TwistMap(E, D);
      q0 := phi(p);
print "q0 = ", q0;
      if Max(Abs(Imaginary(q0[1])), Abs(Imaginary(q0[2]))) gt 10^(-2) then
         q := TracePoint(q0);
      else
         q := q0;
      end if;
print "q = ", q;
      t, Q := Recognize(q);      
   end if;
   if not t then
      print "WARNING: Couldn't correctly recognize point for D=", D;
   end if;

   return EE!p, Q;
end intrinsic;

intrinsic ApproximateRational(x::FldPrElt) -> FldRatElt
{Try to approximate x using continued fractions.}
   vprint Heegner: "Approximating x = ", x;
   cf := ContinuedFraction(Real(x));
   vprint Heegner: "Continued fraction = ", cf;
   for i in [2..#cf] do
      if cf[i] gt 1000 then
         v := Convergents([cf[j] : j in [1..i-1]]);
         return v[1,1]/v[2,1];
      end if;
   end for;
   v := Convergents(cf);
   p := v[1,1]/v[2,1];
   return p;
end intrinsic

intrinsic Recognize(P::PtEll) ->  BoolElt, PtEll
{}
   EC := Scheme(Parent(P));
   E := EllipticCurve([Round(Real(a)) : a in aInvariants(EC)]);
   if P eq EC!0 then
      return true, E!0;
   end if;
   x := ApproximateRational(P[1]);
   R<z> := PolynomialRing(Rationals());
   roots := Roots(Evaluate(DefiningPolynomial(E),[x,z,1]));
   if #roots eq 0 then
      vprint Heegner : "Unable to correctly identify rational point!";
      print "P = ", P;
      return false, x;
   end if;
   if #roots eq 1 then
      return true, E![x,roots[1][1]];
   end if;
   if Abs(roots[1][1]-P[2]) lt Abs(roots[2][1]-P[2]) then
      return true, E![x,roots[1][1]];
   else
      return true, E![x,roots[2][1]];
   end if;
end intrinsic;

intrinsic TwistMap(E::CrvEll, D::RngIntElt) -> Map
{The map phi: E_K --> F_K, where F is the quadratic twist of E
by K = Q(sqrt(D)), and the model for F is Z-minimal.  The second
argument is the base extension of phi to the complex numbers.}
   F := MinimalModel(QuadraticTwist(E, D));
   K := QuadraticField(D);
   C := ComplexField();   
   E_K := BaseExtend(E,K); 
   F_K := BaseExtend(F,K);
   phi_K := Isomorphism(E_K, F_K);
   E_C := EllipticCurve([CoerceIntoC(a) : a in aInvariants(E)]);
   F_C := EllipticCurve([CoerceIntoC(a) : a in aInvariants(F)]);
   r,s,t,u := Explode([CoerceIntoC(a) : a in IsomorphismData(phi_K)]);
   function doit(P, FF) 
      Q := [u^2*P[1]+r, u^3*P[2]+s*u^2*P[1]+t]; 
      if not IsPoint(FF,Q) then
         print "*WARNING*: Possible precision problem!";
         return FF!0;
      end if;
      return FF!Q;
   end function;
   phi_C := hom< E_C -> F_C | P :-> doit(P, F_C)> ;
   return phi_C, [r,s,t,u];
end intrinsic;

intrinsic TracePoint(P::PtEll) -> PtEll
{P + Pbar}
   if P[3] eq 0 then
      return P;
   end if;
   E := Scheme(Parent(P));
   Pbar := E![ComplexConjugate(P[1]), ComplexConjugate(P[2])];
   return P + Pbar;
end intrinsic;


intrinsic Rank0ShaBound(E::CrvEll) 
{Compute an upper bound on Sha for the rank 0 elliptic curve E.}

   A := [-D : D in [1..200] | HeegnerHypothesis(E,-D)]; 
   print "First few Heegner primes: ", A;
   for D in [A[1]] do 
      F := MinimalModel(QuadraticTwist(E,D)); 
      time P := HeegnerPoint(E,D : wp_prec := 100);  
      phi := TwistMap(E,D);
      T := TracePoint(phi(P));
print "T = ", T;
      t, TQ := Recognize(T);  
print "TQ = ", TQ;
      if not t then
         print "Could not recognize T!";
      else
         h := Height(TQ);
         R := Regulator(F);
         s := Sqrt(h/R);
         print "s = ", s;
         a := ApproximateRational(s);
         print "a = ", a;
      end if;
   end for;

end intrinsic;
