CoCalc Public Fileswww / papers / anti-cyclotomic_height_pairing / anti-cyclotomic.mOpen with one click!
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
14
declare verbose ac_height, 2;
15
16
intrinsic 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[2]) and #Factorization(ideal<O|fac[1]>) eq 1 then
45
return false;
46
end if;
47
end for;
48
49
50
// condition (d)
51
E_D := QuadraticTwist(E,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;
57
end intrinsic;
58
59
60
intrinsic 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)]);
73
end intrinsic;
74
75
76
intrinsic 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;
87
end intrinsic;
88
89
90
intrinsic QuadraticIdeal(I::RngOrdIdl, phi::Map) -> RngQuadIdl
91
{The ideal of the appropriate quadratic field defined by the
92
ideal I of a quadratic number field. This function exists only
93
because quadratic fields and number fields are distinct types
94
in MAGMA. Phi is a map from FractionField(Order(I)) to a
95
quadratic field.}
96
"Quad Ideal";
97
vprint ac_height : "QuadraticIdeal()";
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;
105
Parent(gens[1]);
106
I := ideal<O | gens>;
107
"I=",I;
108
return I;
109
end intrinsic;
110
111
112
intrinsic 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[1]+rep[2]*pi)^((p-1)*u);
140
pibar_alpha_v := (rep[1]+rep[2]*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 );
163
end intrinsic;
164
165
intrinsic 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);
177
1;
178
O := MaximalOrder(L);
179
I := ideal<O|1>;
180
xI := x*I;
181
2;
182
time IxI := I meet xI;
183
B := [L!z : z in Basis(IxI)];
184
a,b := Explode(B);
185
3;
186
J := a/x*O + b/x*O;
187
I := ideal<O|a,b>;
188
4;
189
return I, J;
190
191
end intrinsic;
192
193
194
intrinsic 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)[2]);
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[1][1],0);
213
pibar := -Evaluate(fac[2][1],0);
214
215
if P[3] eq 0 then
216
return pAdicField(p)!0;
217
end if;
218
x := P[1]/P[3];
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
230
end intrinsic;
231
232
233
intrinsic 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");
245
end intrinsic;
246
247
248
intrinsic PointsOverK(P::[PtEll], D::RngIntElt, E::CrvEll) -> PtEll
249
{[P] is on Weierstrass model for twist of E by D.}
250
x := [p[1]/p[3] : p in P];
251
y := [p[2]/p[3] : 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[1]/q[3] : q in Q];
261
Qy := [q[2]/q[3] : 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;
269
end intrinsic;
270
271
272
intrinsic PointOverK(P::PtEll, D::RngIntElt, E::CrvEll) -> PtEll
273
{P is on Weierstrass model for twist of E by D.}
274
x := P[1]/P[3];
275
y := P[2]/P[3];
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[1]/Q[3];
285
Qy := Q[2]/Q[3];
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;
293
end intrinsic;
294
295
intrinsic 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);
322
end intrinsic;
323
324
intrinsic 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;
332
end intrinsic;
333
334
335
intrinsic 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;
392
end intrinsic;
393
394
395
396
intrinsic 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
404
return EllipticCurve(CremonaDatabase(),curves[n]), curves[n];
405
end intrinsic;
406
407
408
intrinsic 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;
428
end intrinsic;
429
430
431
intrinsic EquidisTest(curve::MonStgElt, p::RngIntElt, n::RngIntElt, Dstart::RngIntElt, Dstop::RngIntElt)
432
{}
433
// curve := "389A";
434
E := EllipticCurve(CremonaDatabase(),curve);
435
D := FindAcceptableK(E,p,Dstop,Dstart);
436
Qp<q> := pAdicField(p);
437
Qp`SeriesPrinting := 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
476
end intrinsic;
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
/* example usage
494
495
Attach("anti-cyclotomic.m");
496
E := EC("389A");
497
SetVerbose("ac_height",2);
498
p := 3;
499
500
Dp := FindAcceptableK(E,p,-35,-1);
501
P := [* *];
502
for d in Dp do
503
Append(~P, PointFromTwist(E,d));
504
end for;
505
EE := Parent(P[2]);
506
G,f := MordellWeilGroup(E);
507
P0 := EE!Eltseq(f(G.1));
508
P1 := EE!Eltseq(f(G.2));
509
H_2(P0+P[2],p);
510
H_2(P1+P[2],p);
511
*/
512
513