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