Ticket #4849: heegner.py

File heegner.py, 12.1 KB (added by mabshoff, 12 years ago)

This is mostly Magma code

Line 
1"""nodoctest
2Heegner points
3"""
4
5#*****************************************************************************
6#       Copyright (C) 2005 William Stein <wstein@gmail.com>
7#
8#  Distributed under the terms of the GNU General Public License (GPL)
9#
10#    This code is distributed in the hope that it will be useful,
11#    but WITHOUT ANY WARRANTY; without even the implied warranty of
12#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13#    General Public License for more details.
14#
15#  The full text of the GPL is available at:
16#
17#                  http://www.gnu.org/licenses/
18#*****************************************************************************
19
20import ellcurve
21
22"""
23
24
25
26Heegner points package.
27> E :=
28
29> E := EllipticCurve([0,0,1,-1,0]);   // 37A
30> A := [-D : D in [1..100] | HeegnerHypothesis(E,-D)]; A;
31[ -3, -4, -7, -11, -40, -47, -67, -71, -83, -84, -95 ]
32P := [* *]; for D in A do time Y := HeegnerPoint(E,D); Append(~P,Y); end for;
33Time: 0.040
34Time: 0.030
35Time: 0.040
36Time: 0.050
37Time: 0.090
38Time: 0.180
39Time: 0.050
40Time: 0.250
41Time: 0.120
42Time: 0.140
43WARNING: Precision overflow or point is 0.  Use L-series to decide or try
44AssertAttribute(FldPr,"Precision",n) or tweaking parameters.
45Time: 0.270
46
47> for i in [1..10] do Recognize(P[i]); end for;
48(-1 : 0 : 1)
49(1 : -1 : 1)
50(0 : 0 : 1)
51(0 : -1 : 1)
52(1 : -1 : 1)
53(0 : 0 : 1)
54(6 : -15 : 1)
55(0 : -1 : 1)
56(0 : 0 : 1)
57(0 : 0 : 1)
58
59
60* BUG!?
61
62> E := EllipticCurve([ 0, 1, 1, -26, -61 ]);
63> time P := HeegnerPoint(E,-107); P;
64> (3.665647323987715053919669397 + 6.116738926164336860764880743 E-25*i : -0.49999999999999999999999867820523179988 + 9.6625393082244036119397240285848446462*i : 1)
65> IdentifyPoint(P);
66> false 6938055/1892723
67// This time it doesn't identify it as a rational point.  In fact,
68// the coordinates of P generate the Hilbert class field of Q(sqrt(-107)),
69// which is really weird because we just took a trace to get down
70// from the Hilbert class field.  What is going on! 
71// Note: Pete Green's Heegner point package gives the same wrong answer.
72
73"""
74
75"""
76intrinsic tau_to_z(E::CrvEll, tau::FldQuadElt, qexp_prec::RngIntElt) -> FldPrElt
77{}
78   C<i> := ComplexField();
79   tau := CoerceIntoC(tau);
80   if Imaginary(tau) lt 0 then
81      tau := ComplexConjugate(tau);
82   end if;
83   Pi := Pi(C);
84   f,an_over_n := qExpansion(E,qexp_prec);
85   e := Exp(2*Pi*i*tau);
86   V := VectorSpace(C,qexp_prec-1);
87   z := InnerProduct(V!an_over_n, V![e^n : n in [1..qexp_prec-1]]);
88   return z;
89end intrinsic;
90
91intrinsic ideal_to_tau(a::RngQuadIdl, NN::RngQuadIdl, K::FldQuad) -> FldQuadElt
92{Given a nonzero integral fractional ideal a of O_K and an integral
93ideal NN such that O_K/NN = Z/NZ, compute tau in K such that
94(E(tau),C(tau)) = (C/a^(-1), N^(-1)*a^(-1)/a^(-1)).}
95   // reduce ideal: gives equivalent point, since just
96   // multiplies a be a scalar.
97    a := Ideal(Reduction(QuadraticForm(a)));
98
99   // compute I1 and I2 in point (I1, I2/I1)
100   I1 := a^(-1);
101   I2 := I1*NN^(-1);
102
103   // Find a basis for I2 of the form omega_1, omega_2/N for
104   // some basis om1 and om2 of I1.   Such a basis
105   // exists because I2/I1 = Z/NZ.  We compute it using
106   // Smith Normal Form.
107   V := VectorSpace(RationalField(),2);
108   W := VectorSpaceWithBasis([V | Eltseq(b) : b in Basis(I2)]);
109   B1 := Basis(I1);
110   A := MatrixAlgebra(Integers(),2)!&cat[Coordinates(W,V!Eltseq(b)) : b in B1];
111
112   // The smith form S of A, together with unimodular
113   // matrices P and Q such that P * A * Q = S
114   S, P, Q := SmithForm(A);
115   // Apply P to the basis B1 to get om1, om2.
116   om1 := P[1,1]*B1[1] + P[1,2]*B1[2];
117   om2 := P[2,1]*B1[1] + P[2,2]*B1[2];
118
119   // Finally, the point z is 1/om2;
120   z := K!(1/om2);
121   return z;
122end intrinsic;
123
124intrinsic ideal_to_tau(f::QuadBinElt) -> FldQuadElt
125{}
126    D := Discriminant(f);
127    K := QuadraticField(D);
128    sqD := K.1;
129    if D mod 4 eq 0 then
130       sqD := 2*sqD;
131    end if;
132    z := (-f[2] + sqD)/(2*f[1]);
133    return z;
134end intrinsic;
135
136intrinsic levelNforms(N::RngIntElt, D::RngIntElt) -> SeqEnum
137{As in Step 11 of Stephens's paper.}
138    beta := Integers()!Sqrt(Integers(4*N)!D);
139    Q := QuadraticForms(D);
140    G,m := ClassGroup(Q);
141    F := Sort([Reduction(m(g)) : g in G]);
142    levelforms := [Q!1 : i in [1..#F]];
143    known     := [false : i in [1..#F]];
144    found := 0;
145    A := 1;
146    N4 := 4*N;
147    alpha := 0;
148    while found lt #F do
149       B := beta + N4*alpha;
150       C := (B^2 - D)/(N4*A);
151       if Denominator(C) eq 1 then
152          test := Q![Integers()|A*N,B,C];
153          i := Index(F, Reduction(test));
154          if not known[i] then
155             found := found + 1;
156             known[i] := true;
157             levelforms[i] := test;
158          end if;
159       end if;
160       alpha := alpha + 1;
161       if alpha gt A+2 then
162          alpha := 0;
163          A := A + 1;
164       end if;
165    end while;
166
167    return levelforms;
168end intrinsic;
169
170intrinsic HeegnerPointUseIdeals(E::CrvEll, D::RngIntElt,
171                  q_prec::RngIntElt, wp_prec::RngIntElt) -> SeqEnum, CrvEll
172{}
173   require IsFundamentalDiscriminant(D) : "Argument 2 must be a fundamental discriminant.";
174   require D lt 0 : "Argument 2 must be negative.";
175   N := Conductor(E);
176   require GCD(D, 4*N) eq 1 : "D and 4*N must be coprime";
177   require Type(BaseRing(E)) eq FldRat : "Argument 1 must be defined over the rational field.";
178
179   K := QuadraticField(D);
180   OK := MaximalOrder(K);
181   // Find an ideal NN with O/NN = Z/NZ.
182   NN := 1*OK;
183   for p in Factorization(N) do
184      Q := Factorization(p[1]*OK)[1][1];
185      if Degree(Q) ne 1 then
186         require false : "Each prime dividing argument 3 must split completely in Q(sqrt(D)).";
187      end if;
188      NN := NN * Q^p[2];
189   end for;
190   G, m := ClassGroup(OK);
191   vprint Heegner :  "Class number =", #G;
192   taus := [ideal_to_tau(m(g), NN, K) : g in G];
193   vprint Heegner :  "taus =", taus;
194   // taus2 := [ideal_to_tau_2(m(g), NN, K) : g in G];
195   // vprint Heegner :  "taus2 =", taus2;
196   z := &+[tau_to_z(E, tau, q_prec) : tau in taus];
197   vprint Heegner :  "z    =", z;
198   // z2 := &+[tau_to_z(E, tau, q_prec) : tau in taus2];
199   // vprint Heegner :  "z2    =", z2;
200   p := z_to_point(E, z, wp_prec) ;
201   vprint Heegner :  "p = ", p;
202   // p2 := z_to_point(E, z2, wp_prec) ;
203   // vprint Heegner :  "p2 = ", p2;
204   return p;
205end intrinsic;
206
207intrinsic HeegnerHypothesis(E::CrvEll, D::RngIntElt) -> BoolElt
208{True if D satisfies the Heegner Hypotheses, so there is a Heegner
209point attached to discriminant D.}
210   if not IsFundamentalDiscriminant(D) or D ge 0 then
211      return false;
212   end if;
213   N := Conductor(E);   
214   return IsSquare(Integers(4*N)!D) and GCD(N,D) eq 1;
215end intrinsic;
216
217intrinsic RootNumber(E::CrvEll) -> RngIntElt
218{The sign in the functional equation for E, so the analytic rank of E is even
219if and only if the root number is +1.}
220   require Type(BaseRing(E)) eq FldRat : "Argument 1 must be defined over the rational field.";
221   eps := -1;
222/*
223   for p in Factorization(Conductor(E)) do
224   end for;
225*/
226end intrinsic;
227
228intrinsic MyHeegnerPoint(E::CrvEll, D::RngIntElt :
229                        q_prec:= -1, wp_prec := 40, t := 6) -> PtEll, PtEll
230{The Heegner point on E corresponding to D (over the complex numbers), and the a corresponding Q-rational point on
231 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.}
232   require E eq MinimalModel(E) : "Argument 1 must be minimal.";
233   require HeegnerHypothesis(E,D) : "The Heegner Hypotheses must hold.";
234   EE := EllipticCurve([ComplexField()|a : a in aInvariants(E)]);
235   tiny := 10^(-6);
236
237   /*
238   if Abs(LOne(E)) lt tiny then  // E has rank > 0
239      if Abs(LPrimeOne(E)) lt tiny or Abs(LOneChi(E,D)) lt tiny then
240         return EE!0, E!0;
241      end if;
242   else      // E has rank 0, so need twist to have rank at most 1
243      if LPrimeOneChi(E,D) lt tiny then
244         return EE!0, E!0;
245      end if;
246   end if;   
247   */
248
249   N := Conductor(E);
250   require Type(BaseRing(E)) eq FldRat : "Argument 1 must be defined over the rational field.";
251
252   if q_prec le 0 then
253      q_prec := 2*N*Ceiling(Log(-D)/Log(5)) + 100;
254      vprint Heegner : "Using q-expansion to precision", q_prec;
255   end if;
256
257   forms := levelNforms(N, D);
258   taus := [ideal_to_tau(f) : f in forms];
259   vprint Heegner :  "h_K   =", #taus;
260   vprint Heegner :  "tau's =", taus;
261   z := &+[tau_to_z(E, tau, q_prec) : tau in taus];
262   vprint Heegner :  "z     =", z;
263 p := z_to_point(E, z, wp_prec : t := t) ;
264   vprint Heegner :  "p     =", p;
265   if not IsPoint(EE,p) then
266      error "*WARNING*: Precision Overflow!";
267   end if;
268   p := EE!p;
269
270   if Abs(LOne(E)) lt tiny then  // E has rank >=1
271      if Max(Abs(Imaginary(p[1])), Abs(Imaginary(p[2]))) gt 10^(-2) then
272         q := TracePoint(p);
273      else
274         q := p;
275      end if;
276      t, Q := Recognize(q);
277   else   // E has rank 0
278      phi := TwistMap(E, D);
279      q0 := phi(p);
280print "q0 = ", q0;
281      if Max(Abs(Imaginary(q0[1])), Abs(Imaginary(q0[2]))) gt 10^(-2) then
282         q := TracePoint(q0);
283      else
284         q := q0;
285      end if;
286print "q = ", q;
287      t, Q := Recognize(q);     
288   end if;
289   if not t then
290      print "WARNING: Couldn't correctly recognize point for D=", D;
291   end if;
292
293   return EE!p, Q;
294end intrinsic;
295
296intrinsic ApproximateRational(x::FldPrElt) -> FldRatElt
297{Try to approximate x using continued fractions.}
298   vprint Heegner: "Approximating x = ", x;
299   cf := ContinuedFraction(Real(x));
300   vprint Heegner: "Continued fraction = ", cf;
301   for i in [2..#cf] do
302      if cf[i] gt 1000 then
303         v := Convergents([cf[j] : j in [1..i-1]]);
304         return v[1,1]/v[2,1];
305      end if;
306   end for;
307   v := Convergents(cf);
308   p := v[1,1]/v[2,1];
309   return p;
310end intrinsic
311
312intrinsic Recognize(P::PtEll) ->  BoolElt, PtEll
313{}
314   EC := Scheme(Parent(P));
315   E := EllipticCurve([Round(Real(a)) : a in aInvariants(EC)]);
316   if P eq EC!0 then
317      return true, E!0;
318   end if;
319   x := ApproximateRational(P[1]);
320   R<z> := PolynomialRing(Rationals());
321   roots := Roots(Evaluate(DefiningPolynomial(E),[x,z,1]));
322   if #roots eq 0 then
323      vprint Heegner : "Unable to correctly identify rational point!";
324      print "P = ", P;
325      return false, x;
326   end if;
327   if #roots eq 1 then
328      return true, E![x,roots[1][1]];
329   end if;
330   if Abs(roots[1][1]-P[2]) lt Abs(roots[2][1]-P[2]) then
331      return true, E![x,roots[1][1]];
332   else
333      return true, E![x,roots[2][1]];
334   end if;
335end intrinsic;
336
337intrinsic TwistMap(E::CrvEll, D::RngIntElt) -> Map
338{The map phi: E_K --> F_K, where F is the quadratic twist of E
339by K = Q(sqrt(D)), and the model for F is Z-minimal.  The second
340argument is the base extension of phi to the complex numbers.}
341   F := MinimalModel(QuadraticTwist(E, D));
342   K := QuadraticField(D);
343   C := ComplexField();   
344   E_K := BaseExtend(E,K);
345   F_K := BaseExtend(F,K);
346   phi_K := Isomorphism(E_K, F_K);
347   E_C := EllipticCurve([CoerceIntoC(a) : a in aInvariants(E)]);
348   F_C := EllipticCurve([CoerceIntoC(a) : a in aInvariants(F)]);
349   r,s,t,u := Explode([CoerceIntoC(a) : a in IsomorphismData(phi_K)]);
350   function doit(P, FF)
351      Q := [u^2*P[1]+r, u^3*P[2]+s*u^2*P[1]+t];
352      if not IsPoint(FF,Q) then
353         print "*WARNING*: Possible precision problem!";
354         return FF!0;
355      end if;
356      return FF!Q;
357   end function;
358   phi_C := hom< E_C -> F_C | P :-> doit(P, F_C)> ;
359   return phi_C, [r,s,t,u];
360end intrinsic;
361
362intrinsic TracePoint(P::PtEll) -> PtEll
363{P + Pbar}
364   if P[3] eq 0 then
365      return P;
366   end if;
367   E := Scheme(Parent(P));
368   Pbar := E![ComplexConjugate(P[1]), ComplexConjugate(P[2])];
369   return P + Pbar;
370end intrinsic;
371
372
373intrinsic Rank0ShaBound(E::CrvEll)
374{Compute an upper bound on Sha for the rank 0 elliptic curve E.}
375
376   A := [-D : D in [1..200] | HeegnerHypothesis(E,-D)];
377   print "First few Heegner primes: ", A;
378   for D in [A[1]] do
379      F := MinimalModel(QuadraticTwist(E,D));
380      time P := HeegnerPoint(E,D : wp_prec := 100); 
381      phi := TwistMap(E,D);
382      T := TracePoint(phi(P));
383print "T = ", T;
384      t, TQ := Recognize(T); 
385print "TQ = ", TQ;
386      if not t then
387         print "Could not recognize T!";
388      else
389         h := Height(TQ);
390         R := Regulator(F);
391         s := Sqrt(h/R);
392         print "s = ", s;
393         a := ApproximateRational(s);
394         print "a = ", a;
395      end if;
396   end for;
397
398end intrinsic;
399"""