Author: William A. Stein
Compute Environment: Ubuntu 18.04 (Deprecated)
1/*******************************************************************
2*
3*   anti-cyclotomic.m   -- test Mazur's "sign conjecture"
4*
5*   William A. Stein
6*
7*   02/28/03, 04/01 (initial creation date).
8*
9*   The file anti-cyclotomic_height.tex, which Mazur wrote,
10*   describes the algorithm that is implemented below.
11*
12********************************************************************/
13
14declare verbose ac_height, 2;
15
16intrinsic IsAcceptableK(E::CrvEll, p::RngIntElt, D::RngIntElt) -> BoolElt
17//   {True if and only if K=Q(sqrt(D)) satisfies conditions (a)--(d).}
18   {True if rank(E) is odd, which for all practical purposes in our cases means 1.}
19
20   vprint ac_height : "IsAcceptableK()";
21
22   if not IsFundamental(D) then
23      return false;
24   end if;
25
26   // condition (a)
27   N := Conductor(E);
28   if GCD(D,N*p) ne 1 then
29      return false;
30   end if;
31
32   // condition (b)
33   R<x> := PolynomialRing(RationalField());
34   K := NumberField(x^2-D);
35   O := MaximalOrder(K);
36   if #Factorization(ideal<O|p>) eq 1 then
37      return false;
38   end if;
39
40
41   // condition (c)   // equiv to condition (d)?
42   // all primes dividing N split
43   for fac in Factorization(N) do
44      if IsOdd(fac) and #Factorization(ideal<O|fac>) eq 1 then
45         return false;
46      end if;
47   end for;
48
49
50   // condition (d)
52   if Evaluate(KroneckerCharacter(D),-N) eq 1 then    // there's no chance to have rank 1
53      return false;
54   end if;
55
56   return true;
57end intrinsic;
58
59
60intrinsic FindAcceptableK(E::CrvEll, p::RngIntElt, Dmin::RngIntElt, Dmax::RngIntElt) -> SeqEnum
61   {The fields K whose discriminant lie between Dmin and Dmax and satisfy conditions (a)--(d).}
62
63   vprint ac_height : "FindAcceptableK()";
64
65   require IsPrime(p) and p gt 0 : "Argument 2 must be a prime number.";
66   require Dmin lt 0: "Argument 3 must be negative.";
67   require Dmax lt 0: "Argument 4 must be negative.";
68   require Dmin le Dmax: "Argument 3 must be less than or equal to argument 4.";
69
70   N := Conductor(E);
71   vprint ac_height : "Finding acceptable K for E of conductor", N;
72   return Reverse([D : D in [Dmin..Dmax] | IsAcceptableK(E,p,D)]);
73end intrinsic;
74
75
76intrinsic FindAcceptableK(E::CrvEll, p::RngIntElt) -> RngIntElt
77{}
78   vprint ac_height : "FindAcceptableK()";
79
80   D := -1;
81   while true do
82      if IsAcceptableK(E,p,D) then
83         return D;
84      end if;
85      D := D - 1;
86   end while;
87end intrinsic;
88
89
91{The ideal of the appropriate quadratic field defined by the
92ideal I of a quadratic number field.  This function exists only
93because quadratic fields and number fields are distinct types
94in MAGMA.   Phi is a map from FractionField(Order(I)) to a
98
99   K := FieldOfFractions(Order(I));
100   L := Codomain(phi);
101   O := MaximalOrder(L);
102"O=",O;
103   gens := [O!phi(K!g) : g in Generators(I)];
104"gens=",gens;
105Parent(gens);
106   I := ideal<O | gens>;
107"I=",I;
108   return I;
109end intrinsic;
110
111
112intrinsic rhotilde(I::RngOrdIdl , pi::RngLocElt, pibar::RngLocElt) -> FldLocElt
113   {I is an ideal of O_K that is coprime to p. Each of pi and pibar define
114   a map from O_K to Z_p; pi is the image of O.1, and similarly for pibar.
115   (Note that pi and pibar are computed by finding the roots of the charpoly
116   of O.1 over Z_p.)}
117
118   vprint ac_height : "rhotilde(I,", pi, ",", pibar,")";
119
120   p := Prime(Parent(pi));
121   time n := Numerator(Norm(I)/1);
122   if n mod p eq 0 then
123      vprint ac_height : "Bad point.";
124      return false;
125   end if;
126   O := Order(I);
127   u := #UnitGroup(O);
128   vprint ac_height: "Computing class number h";
129   h := ClassNumber(O);
130   vprint ac_height: "Class number h =", h;
131   v := (p-1)*u*h;
132   vprint ac_height : "Raising ideal to the power h";
133   Ih  := I^h;
134   vprint ac_height : "Finding generator for the principal ideal I^h";
135   _, alpha_h := IsPrincipal(Ih);
136   vprint ac_height : "I^h is principal.";
137   rep := Eltseq(alpha_h);
138   vprint ac_height : "Computing embeddings of canonical power of generator into Qp";
139   pi_alpha_v    := (rep+rep*pi)^((p-1)*u);
140   pibar_alpha_v := (rep+rep*pibar)^((p-1)*u);
141   //vprint ac_height : "pi_alpha_v =", pi_alpha_v;
142   //vprint ac_height : "pibar_alpha_v =", pibar_alpha_v;
143
144   //vprint ac_height : "pi_alpha_v-1 = ", pi_alpha_v-1;
145   val1 := Valuation(pi_alpha_v-1);
146   // vprint ac_height : "Valuation(pi_alpha_v-1) = ", val1;
147   //vprint ac_height : "pibar_alpha_v-1 = ", pibar_alpha_v-1;
148   val2 := Valuation(pibar_alpha_v-1);
149   //vprint ac_height : "Valuation(pibar_alpha_v-1) = ", val2;
150
151   if val1*val2 eq 0 then
152      vprint ac_height : "Bad point.";
153      return false;
154   end if;
155
156   vprint ac_height : "Valuations ok -- now computing p-adic logarithms.";
157   logpi := Log(pi_alpha_v);
158   //vprint ac_height : "lambda(pi_alpha_v) =", logpi;
159   logpibar := Log(pibar_alpha_v);
160   //vprint ac_height : "lambda(pibar_alpha_v) =", logpibar;
161
162   return 1/v * ( logpi - logpibar );
163end intrinsic;
164
165intrinsic Num_and_Denom_Ideals(x::FldNumElt) -> RngOrdIdl
166   {I and J, where I is the numerator ideal of x and J is the denominator ideal of x.}
167
168   // We compute the denominator ideal as follows.  Compute
169   // the intersection of the image of multiplication by x
170   // and O inside K.  Let a and b be generators for this intersection,
171   // as a Z-module.  Then the denominator ideal is generated by
172   // a/x and b/x and the numerator ideal is generated by a and b.
173
174   vprint ac_height : "Num_and_Denom_Ideals()";
175
176   L   := Parent(x);
1771;
178   O   := MaximalOrder(L);
179   I   := ideal<O|1>;
180   xI  := x*I;
1812;
182   time IxI := I meet xI;
183   B   := [L!z : z in Basis(IxI)];
184   a,b := Explode(B);
1853;
186   J := a/x*O + b/x*O;
187   I := ideal<O|a,b>;
1884;
189   return I, J;
190
191end intrinsic;
192
193
194intrinsic H_2(P::PtEll, p::RngIntElt) -> FldLocElt
195   {The function H_2 defined in Mazur's note.  This is rhotilde(J*x(P))/2, where
196   J is the "ideal denominator" of the principal ideal (x(P)) generated by
197   the x-coordinate of P.}
198
199   vprint ac_height : "H_2()";
200
201   require IsPrime(p) and p gt 0 : "Argument 2 must be a prime number.";
202
203   E := Curve(Parent(P));
204   K := BaseRing(E);
205   require Type(K) eq FldNum and Degree(K) eq 2 and Discriminant(K) lt 0 :
206         "The base field of the parent of argument 1 must be a quadratic imaginary number field.";
207
208   O := MaximalOrder(K);
209   f := PolynomialRing(pAdicRing(p : Precision := 40))!MinimalPolynomial(Basis(O));
210   fac := Factorization(f);
211   require #fac eq 2 : "Argument 2 must split in the base field of the parent of argument 1.";
212   pi := -Evaluate(fac,0);
213   pibar := -Evaluate(fac,0);
214
215   if P eq 0 then
217   end if;
218   x := P/P;
219   if Numerator(Norm(x)) mod p eq 0 then
220      error "Can't compute height of point.";
221   end if;
222
223   I, J := Num_and_Denom_Ideals(x);
224   r := rhotilde(I, pi, pibar);
225   if Type(r) eq BoolElt then
226      return r;
227   end if;
228   return r/2;
229
230end intrinsic;
231
232
233intrinsic Pairing(P::PtEll, Q::PtEll, p::RngIntElt, n::RngIntElt) -> FldLocElt
234   {Let h(x) = H_2(p^n*x)/p^(2n).  Then this returns
235   h(P+Q) + h(P-Q) - 2*h(P) - 2*h(Q),
236   which is an approximation for the height pairing of P and Q.}
237
238   vprint ac_height : "Pairing";
239
240   function h(x,p, comment)
241      vprint ac_height : "Computing height of ", comment, " using n = ", n;
242      return H_2(p^n*x,p)/p^(2*n);
243   end function;
244   return h(P+Q,p, "P+Q") + h(P-Q,p,"P-Q") - 2*h(P,p,"P") - 2*h(Q,p,"Q");
245end intrinsic;
246
247
248intrinsic PointsOverK(P::[PtEll], D::RngIntElt, E::CrvEll) -> PtEll
249{[P] is on Weierstrass model for twist of E by D.}
250   x := [p/p : p in P];
251   y := [p/p : p in P];
252
253   R<z> := PolynomialRing(Rationals());
254   K := NumberField(z^2-D);
255   sqrtD := K.1;
256   Ew,_,psi := WeierstrassModel(E);
257   Ew := BaseExtend(Ew,K);
258   Pw := [[x[i]/D, (y[i]/D^2)*sqrtD] : i in [1..#P]];
259   Q := [Ew!p : p in Pw];
260   Qx := [q/q : q in Q];
261   Qy := [q/q : q in Q];
262   r,s,t,u := Explode(IsomorphismData(psi));
263
264   // (x, y) |-> (u^2x + r, u^3 y + su^2 x + t).
265   X := [u^2*Qx[i] + r : i in [1..#P]];
266   Y := [u^3*Qy[i] + s*u^2*Qx[i] + t : i in [1..#P]];
267   E := BaseExtend(E,K);
268   return [E![X[i],Y[i]] : i in [1..#P]], E;
269end intrinsic;
270
271
272intrinsic PointOverK(P::PtEll, D::RngIntElt, E::CrvEll) -> PtEll
273{P is on Weierstrass model for twist of E by D.}
274   x := P/P;
275   y := P/P;
276
277   R<z> := PolynomialRing(Rationals());
278   K := NumberField(z^2-D);
279   sqrtD := K.1;
280   Ew,_,psi := WeierstrassModel(E);
281   Ew := BaseExtend(Ew,K);
282   Pw := [x/D, (y/D^2)*sqrtD];
283   Q := Ew!Pw;
284   Qx := Q/Q;
285   Qy := Q/Q;
286   r,s,t,u := Explode(IsomorphismData(psi));
287
288   // (x, y) |-> (u^2x + r, u^3 y + su^2 x + t).
289   X := u^2*Qx + r;
290   Y := u^3*Qy + s*u^2*Qx + t;
291   E := BaseExtend(E,K);
292   return E![X,Y], E;
293end intrinsic;
294
295intrinsic PointFromTwist(E::CrvEll, D::RngIntElt) -> PtEll
296   {The point on E/Q(sqrt(D)) corresponding to a generator of E_D(Q)/tor,
297   assuming that the latter group has rank 1.}
298
299   vprint ac_height : "PointFromTwist()";
300
301   Ew := WeierstrassModel(E);
302   _,_,_,a,b  := Explode(aInvariants(Ew));
303   E_D := EllipticCurve([D^2*a,D^3*b]);
304   G,f := MordellWeilGroup(E_D : Bound := 2, HeightBound := 0.5);
305   //   require TorsionFreeRank(G) eq 1 : "The twist of argument 1 by argument 2 must have rank 1.";
306   if TorsionFreeRank(G) gt 1 then
307      print "WARNING: rank is > 1 for E = ", E, " D = ", D;
308   end if;
309   if TorsionFreeRank(G) lt 1 then
310      return false;
311      /*
312      print "WARNING: rank of E_D computed is < 1.  We have E = ", aInvariants(E), " D = ", D;
313      print "WARNING: rank computed is < 1 for E_D = ", aInvariants(E_D);
314      G,f := MordellWeilGroup(E_D);  // this might take longer...!
315      print "Computed generators and go = ", G;
316      */
317   end if;
318
319   // P is the generator for E_D(Q)/tor
320   P := f(G.Ngens(G));
321   return PointOverK(P, D, E);
322end intrinsic;
323
324intrinsic TestNontrivialityConjecture(E::CrvEll, p::RngIntElt,
325                                      n::RngIntElt) -> List, RngIntElt
326{}
327   require IsPrime(p) and p gt 0 : "Argument 2 must be a prime number.";
328   require Conductor(E) mod p ne 0 : "Argument 1 must have good reduction at argument 2.";
329
330   D := FindAcceptableK(E,p);
331   return TestNontrivialityConjecture(E,p,n,[D]), D;
332end intrinsic;
333
334
335intrinsic TestNontrivialityConjecture(E::CrvEll, p::RngIntElt, n::RngIntElt, Dlist::SeqEnum)
336       -> List
337{Let E be a rank two elliptic curve and let P1, P2 be a basis for E(Q).
338 Let P = P1 + P0, where P0 is a point of infinite order in E_D(Q) subset E(Q(sqrt(D))).
339 The value returned is a list of pairs
340
341          <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) *]>
342
343 Some of the H_2 values might equal false, if the point doesn't like in the
344 appropriate subgroup for H_2 to be defined.
345}
346
347   require IsPrime(p) and p gt 0 : "Argument 2 must be a prime number.";
348   require Conductor(E) mod p ne 0 : "Argument 1 must have good reduction at argument 2.";
349
350   G,f := MordellWeilGroup(E);
351   require TorsionFreeRank(G) eq 2 : "Argument 1 must have rank 2.";
352   P1 := f(G.(Ngens(G)-1));
353   P2 := f(G.Ngens(G));
354
355   ans := [* *];
356   for D in Dlist do
357      P0 := PointFromTwist(E,D);
358      if Type(P0) eq BoolElt then
359	 print "Couldn't find a point on the twist quickly enough, so bailing.";
360	 print "E = ", E;
361	 print "D = ", D;
362	 continue;
363      end if;
364      Pa  := P0 + Parent(P0)!Eltseq(P1);
365      Pb  := P0 + Parent(P0)!Eltseq(P2);
366      hts := [* *];
367      for m in [1..n] do
368         H2a := H_2(Pa,p);
369         H2b := H_2(Pb,p);
370
371         if Type(H2a) ne BoolElt then
372            H2a := H2a/p^(2*m);
373         end if;
374         if Type(H2b) ne BoolElt then
375            H2b := H2b/p^(2*m);
376         end if;
377         Append(~hts, <H2a, H2b>);
378
379         vprint ac_height : "When m =", m, " hts = ", hts;
380
381         if m lt n then
382            vprint ac_height : "Multiplying each point by", p;
383            time Pa := p*Pa;
384            time Pb := p*Pb;
385         end if;
386      end for;
387      Append(~ans, <D, hts>);
388      vprint ac_height : "Answer so far = ", ans;
389   end for;
390
391   return ans;
392end intrinsic;
393
394
395
396intrinsic Rank2Curve(n::RngIntElt) -> CrvEll
397{}
398   curves := ["389A","433A","446D", "563A","571B","643A","655A","664A","681C",
399              "707A","709A","718B","794A","817A","916C","944E","997B","997C","1001C",
400              "1028A","1034A","1058C","1070A","1073A","1077A","1088J","1094A","1102A",
401              "1126A","1132A","1137A","1141A","1143C","1147A","1171A"];
402   require n ge 1 and n le #curves : "Argument 1 must be between 1 and", #curves;
403
405end intrinsic;
406
407
408intrinsic SmallestDTest(p::RngIntElt)
409{}
410   out := Open(Sprintf("SmallestDTest_p%o.out",p),"w");
411
412   for n in [1..35] do
413      E, label := Rank2Curve(n);
414      if Conductor(E) mod p eq 0 then
415         continue;
416      end if;
417      f := qEigenform(E,p+1);
418      if Coefficient(f,p) eq 0 then
419	 continue;
420      end if;
421      fprintf out, "\n// ** %o **: %o\n", label, E;
422      Flush(out);
423      ans, D := TestNontrivialityConjecture(E,p,4);
424      fprintf out, "\n// Q(sqrt(%o))\n\n", D;
425      fprintf out, "heights[%o] := %o;\n\n", n, ans;
426      Flush(out);
427   end for;
428end intrinsic;
429
430
431intrinsic EquidisTest(curve::MonStgElt, p::RngIntElt, n::RngIntElt, Dstart::RngIntElt, Dstop::RngIntElt)
432{}
433//   curve := "389A";
435   D := FindAcceptableK(E,p,Dstop,Dstart);
437   QpSeriesPrinting := true;
438
439   vprint ac_height : "Finding mordell Weil group of", E;
440   G,f := MordellWeilGroup(E);
441   out := Open(Sprintf("equidis_test_%o_p%o_n%o_d%o.m", curve, p, n, Dstart),"w");
442   fprintf out, "D = %o\n", D;
443
444   for d in D do
445      t := Cputime();
446      fprintf out, "\nK = Q(sqrt(%o)):\t", d;
447      vprint ac_height : "Finding point on twist by", d;
448      P  := PointFromTwist(E,d);
449      vprint ac_height : "The point is P = ", P;
450      if Type(P) eq BoolElt then
451         continue;
452      end if;
453      EE := Parent(P);
454      P0 := EE!Eltseq(f(G.1));
455      P1 := EE!Eltseq(f(G.2));
456      vprint ac_height : "Computing first height.";
457      hAC := H_2(p^n*(P0+P),p);
458      if Type(hAC) ne BoolElt then
459         hAC := hAC/p^(2*n);
460         vprint ac_height : "Computing second height.";
461         hBC := H_2(p^n*(P1+P),p);
462         if Type(hBC) ne BoolElt then
463            hBC := hBC/p^(2*n);
464         end if;
465      end if;
466      if Type(hAC) eq BoolElt or Type(hBC) eq BoolElt then
467         fprintf out, "(no data)";
468      else
469         r := hAC/hBC;
470         fprintf out, "%o", r;
471      end if;
472      fprintf out, " (time = %o)\n", Cputime(t);
473      Flush(out);
474   end for;
475
476end intrinsic;
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493/* example usage
494
495Attach("anti-cyclotomic.m");
496E := EC("389A");
497SetVerbose("ac_height",2);
498p := 3;
499
500Dp := FindAcceptableK(E,p,-35,-1);
501P := [* *];
502for d in Dp do
503   Append(~P, PointFromTwist(E,d));
504end for;
505EE := Parent(P);
506G,f := MordellWeilGroup(E);
507P0 := EE!Eltseq(f(G.1));
508P1 := EE!Eltseq(f(G.2));
509H_2(P0+P,p);
510H_2(P1+P,p);
511*/
512
513
`