CoCalc Public Fileswww / talks / harvard-talk-2004-12-08 / padic_height.m
Author: William A. Stein
Compute Environment: Ubuntu 18.04 (Deprecated)
1
/*****************************************************************
2
3
4
cyclotomic height pairing.
5
6
(c) 2004, William Stein, was@math.harvard.edu
7
8
****************************************************************/
9
10
/*****************************************************************
11
1. You *must* put a copy of kedlaya.m in the same directory as this file,
12
and attach it with
13
14
> Attach("kedlaya.m");
15
16
The following line imports five functions from the kedlaya.m
17
that is distributed with MAGMA (as package/Geometry/CrvHyp/kedlaya.m)
18
The functions FrobYInv, Convert, ReduceA, and ReduceB implement
19
arithmetic in Monsky-Washnitzer cohomology, following Kedlaya's
20
algorithms. These functions are used below in the intrinsic E2
21
but nowhere else.
22
****************************************************************/
23
24
import "kedlaya.m": myXGCD, FrobYInv, Convert, ReduceA, ReduceB;
25
26
/*****************************************************************
27
2. How do you use this package?
28
29
First try the following calculations and see if they work:
30
31
> Attach("kedlaya.m");
32
33
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
34
> P := E![0,0];
35
> h := height_function(E, 5, 20);
36
> h(P);
37
-201147061754703 + O(5^21)
38
39
Also, compute some regulators:
40
41
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
42
> regulator(E,5,20);
43
138360994885922 + O(5^21)
44
> regulator(E,7,20);
45
4709403600911866 + O(7^20)
46
47
> regulator(E,7,20);
48
22975764581280320 + O(7^20)
49
50
Now you should know enough to be able to make use of the package.
51
However, you should probably read through the source code below, To
52
make this easier, I've put a description of *every* function with a
53
complete example of usage. I've also grouped the functions into
54
logical sections.
55
56
****************************************************************/
57
58
59
/*****************************************************************
60
61
Formal groups of elliptic curves.
62
63
****************************************************************/
64
65
R<t> := LaurentSeriesRing(RationalField());
66
67
intrinsic formal_w(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
68
{
69
Compute and return the formal power series w(t) defined in Chapter IV, Section 1
70
of Silverman's AEC book. This function is useful because Silverman gives
71
expressions in w(t) for other formal power series related to formal groups.
72
73
EXAMPLE:
74
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
75
> formal_w(E,19);
76
t^3 + t^6 - t^7 + 2*t^9 - 4*t^10 + 2*t^11 + 5*t^12 - 15*t^13 +
77
15*t^14 + 9*t^15 - 56*t^16 + 84*t^17 - 14*t^18 + O(t^19)
78
}
79
w := t^3 + O(t^prec);
80
a1, a2, a3, a4, a6 := Explode(aInvariants(E));
81
for n in [4..prec-1] do
82
w := w + t^n*(a1*Coefficient(w,n-1) +
83
a2*Coefficient(w,n-2) +
84
a3*&+[Integers()|Coefficient(w,i)*Coefficient(w,n-i) : i in [3..n-3]] +
85
a4*&+[Integers()|Coefficient(w,i)*Coefficient(w,n-1-i) : i in [3..n-4]] +
86
a6*&+[Integers()|Coefficient(w,i)*Coefficient(w,j)*Coefficient(w,n-i-j) :
87
i in [3..n-3], j in [3..n-3] | n-i-j ge 3]);
88
end for;
89
return w;
90
end intrinsic;
91
92
intrinsic formal_x(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
93
{
94
Compute and return the formal group power series defined by the
95
function x on E.
96
97
EXAMPLES:
98
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
99
> formal_x(E,15);
100
t^-2 - t + t^2 - t^4 + 2*t^5 - t^6 - 2*t^7 + 6*t^8 - 6*t^9 -
101
3*t^10 + 20*t^11 - 30*t^12 + 6*t^13 + 65*t^14 + O(t^15)
102
}
103
return t/formal_w(E,prec+5);
104
end intrinsic;
105
106
intrinsic formal_y(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
107
{
108
Compute and return the formal group power series defined by the
109
function y on E.
110
111
EXAMPLES:
112
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
113
> formal_y(E,15);
114
-t^-3 + 1 - t + t^3 - 2*t^4 + t^5 + 2*t^6 - 6*t^7 + 6*t^8 + 3*t^9
115
- 20*t^10 + 30*t^11 - 6*t^12 - 65*t^13 + 140*t^14 + O(t^15)
116
}
117
return -1/formal_w(E,prec+6);
118
end intrinsic;
119
120
121
intrinsic formal_omega(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
122
{
123
Compute and return the formal differential in terms of t, i.e., return
124
the function f(t) such that f(t)dt is the formal differential omega.
125
126
ALGORITHM:
127
Use that omega = dx/(2*y + a1*x + a3) = f(t) dt.
128
129
EXAMPLES:
130
> formal_omega(E,10);
131
1 + 2*t^3 - 2*t^4 + 6*t^6 - 12*t^7 + 6*t^8 + 20*t^9 - 60*t^10 +
132
60*t^11 + O(t^12)
133
}
134
135
a1, a2, a3, a4, a6 := Explode(aInvariants(E));
136
x := formal_x(E,prec);
137
y := formal_y(E,prec);
138
139
return Derivative(x)/(2*y + a1*x + a3);
140
141
end intrinsic;
142
143
intrinsic formal_log(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
144
{
145
The formal log in terms of t. This is sometimes written z=z(s). It is the
146
integral with respect to t of the formal differential.
147
148
EXAMPLE:
149
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
150
> formal_log(E,10);
151
t + 1/2*t^4 - 2/5*t^5 + 6/7*t^7 - 3/2*t^8 + 2/3*t^9 + 2*t^10 -
152
60/11*t^11 + 5*t^12 + O(t^13)
153
}
154
y := formal_y(E,prec);
155
x := formal_x(E,prec);
156
a1, a2, a3, a4, a6 := Explode(aInvariants(E));
157
return Integral(1/(2*y + a1*x + a3) * Derivative(x));
158
end intrinsic;
159
160
intrinsic formal_inverse(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
161
{
162
Compute the power series in Z[[t]] corresponding to the inverse
163
map in the formal group on E.
164
165
EXAMPLES:
166
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
167
> formal_inverse(E,10);
168
-t - t^4 - 2*t^7 + t^8 - 5*t^10 + 6*t^11 - 2*t^12 + O(t^13)
169
}
170
x := formal_x(E,prec);
171
y := formal_y(E,prec);
172
a1, a2, a3, a4, a6 := Explode(aInvariants(E));
173
return x/(y + a1*x + a3);
174
end intrinsic;
175
176
intrinsic formal_group_law(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
177
{
178
Compute and return the formal group law on E as an element of Z[[t,t2]].
179
The precision must be at least 4.
180
181
ALGORITHM:
182
Use the formula at the end of Section 1 of Chapter 4 of Silverman AEC,
183
on page 114. This is the formula for "z_3" there. Note though that
184
the FORMULA IS *WRONG* in Silverman's book!! The a_1 and a_3 in that
185
formula should have minus signs in front, but they don't.
186
187
188
EXAMPLE:
189
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
190
> formal_group_law(E,4);
191
t + O(t^4) + (1 - 2*t^3 + 2*t^4 + O(t^6))*t2 + (-3*t^2 + 4*t^3 +
192
O(t^5))*t2^2 + (-2*t + 4*t^2 + O(t^4))*t2^3 + O(t2^4)
193
}
194
require prec ge 4: "The precision must be at least 4.";
195
t1 := t;
196
S<t2> := LaurentSeriesRing(R);
197
w := formal_w(E,prec);
198
function tsum(n)
199
return &+[t2^m * t1^(n-m-1) : m in [0..n-1]];
200
end function;
201
lam := &+[tsum(n)*Coefficient(w,n) : n in [3..prec-1]];
202
203
w1 := Evaluate(w,t1+O(t^prec));
204
nu := w1 - lam*(t1+O(t^prec));
205
a1, a2, a3, a4, a6 := Explode(aInvariants(E));
206
t3 := -t1 - t2 + \
207
(-a1*lam - a3*lam^2 - a2*nu - 2*a4*lam*nu - 3*a6*lam^2*nu)/
208
(1 + a2*lam + a4*lam^2 + a6*lam^3) + O(t1^prec) + O(t2^prec);
209
inv := formal_inverse(E, prec);
210
return Evaluate(inv, t3);
211
end intrinsic;
212
213
intrinsic formal_mult(E::CrvEll, n::RngIntElt, prec::RngIntElt) -> RngSerLaurElt
214
{
215
Compute and return the formal power series corresponding to multiplication by n on E,
216
to precision prec.
217
218
EXAMPLE:
219
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
220
> formal_mult(E,2,5);
221
2*t - 7*t^4 + 12*t^5 + 4*t^7 + O(t^8)
222
> formal_mult(E,3,5);
223
3*t - 39*t^4 + 96*t^5 + 234*t^7 + O(t^8)
224
> formal_mult(E,7,5);
225
7*t - 1197*t^4 + 6720*t^5 + 115254*t^7 + O(t^8)
226
}
227
228
if n eq -1 then
229
return formal_inverse(E,prec);
230
end if;
231
232
// The following is less direct, but it works more quickly...
233
L := FunctionField(E);
234
EL := BaseExtend(E,L);
235
Q := EL![L.1,L.2];
236
x,y := Explode(Eltseq(n*Q));
237
a := -x/y;
238
xx := formal_x(E,prec);
239
yy := formal_y(E,prec);
240
phi := hom<Parent(x) -> Parent(xx) | xx,yy>;
241
return phi(a);
242
end intrinsic;
243
244
intrinsic formal_divpoly(E::CrvEll, m::RngIntElt, prec::RngIntElt) -> RngSerElt
245
{
246
Compute and return the formal division "polynomial" f_m for E to
247
precision prec.
248
249
EXAMPLES:
250
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
251
> formal_divpoly(E,2,5);
252
2*t^-3 - 3 + 2*t - 2*t^3 + O(t^4)
253
> formal_divpoly(E,3,5);
254
3*t^-8 - 12*t^-5 + 6*t^-4 + 9*t^-2 + O(t^-1)
255
> formal_divpoly(E,7,5);
256
7*t^-48 - 168*t^-45 - 140*t^-44 + 2750*t^-42 + O(t^-41)
257
}
258
if m eq 1 then
259
return R!1;
260
end if;
261
if m eq -1 then
262
return R!(-1);
263
end if;
264
if m eq 2 then
265
return Sqrt(Evaluate(DivisionPolynomial(E,m),formal_x(E,prec)));
266
end if;
267
assert IsOdd(m);
268
return Evaluate(DivisionPolynomial(E,m),formal_x(E,prec));
269
end intrinsic;
270
271
intrinsic weierstrass_p(E::CrvEll, prec::RngIntElt) -> FldPrElt
272
{
273
Compute and return the Weierstrass p-series for E to precision prec.
274
275
EXAMPLES:
276
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
277
> weierstrass_p(E,9);
278
t^-2 + 1/5*t^2 - 1/28*t^4 + 1/75*t^6 - 3/1540*t^8 + O(t^10)
279
}
280
L<w> := PolynomialRing(Rationals());
281
R<t> := LaurentSeriesRing(L);
282
p := t^(-2);
283
c4, c6 := Explode(cInvariants(E));
284
a := 4; b := -c4/12; c := -c6/216;
285
for n in [1..prec] do
286
ptest := p + w*t^n + O(t^(n+1));
287
ptest_prime := Derivative(ptest);
288
f := ptest_prime^2 - a*ptest^3 - b*ptest - c;
289
x := Coefficient(f,-4+n);
290
coeff := -Coefficient(x,0)/Coefficient(x,1);
291
p := p + coeff*t^n;
292
end for;
293
return LaurentSeriesRing(Rationals())!(p + O(t^(prec+1)));
294
end intrinsic;
295
296
297
intrinsic t_val(Q::PtEll) -> FldRatElt
298
{
299
Given a point (x,y) on an elliptic curve, return t = -x/y. (This function doesn't
300
do much!)
301
302
EXAMPLE:
303
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
304
> P := E![0,0]; P;
305
(0 : 0 : 1)
306
> Q := 5*P;
307
(1 : 0 : 1)
308
> Q := 5*P;
309
> Q;
310
(1/4 : -5/8 : 1)
311
> t_val(Q);
312
2/5
313
}
314
require Q[2] ne 0 : "The denominator y must be nonzero.";
315
return -Q[1]/Q[2];
316
end intrinsic;
317
318
319
/*****************************************************************
320
321
Computing sigma: Integrality algorithm
322
323
****************************************************************/
324
325
intrinsic formal_sigma_in_s2(E::CrvEll, prec::RngIntElt) -> RngSerElt
326
{
327
Compute a formal power series for sigma with coefficient polynomials
328
in the coefficient s_2 of t^3.
329
330
NOTE: This is adapted from Christian Wuthrich's PARI code, which is also the
331
algorithm Mazur-Tate and I found on 2004-06-09.
332
333
EXAMPLE:
334
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
335
> formal_sigma_in_s2(E, 10);
336
t + s2*t^3 + 1/2*t^4 + (1/2*s2^2 - 5/12)*t^5 + 3/2*s2*t^6 +
337
(1/6*s2^3 - 73/60*s2 + 103/120)*t^7 + (5/4*s2^2 - 37/24)*t^8
338
+ (1/24*s2^4 - 121/120*s2^2 + 2791/840*s2 + 1411/2016)*t^9 +
339
(7/12*s2^3 - 691/120*s2 + 481/240)*t^10 + (-7/15*s2^3 +
340
95/28*s2^2 + 379/150*s2 - 85793/15400)*t^11 + (3/16*s2^4 -
341
463/80*s2^2 + 4873/560*s2 + 6977/1344)*t^12 + O(t^13)
342
}
343
344
R<s2> := FieldOfFractions(PolynomialRing(RationalField()));
345
S<z> := LaurentSeriesRing(R);
346
w := S!weierstrass_p(E,prec);
347
w := &+[Coefficient(w,n)*z^n : n in [0..prec-1]] + O(z^prec);
348
sigma_of_z := z*Exp( s2*z^2 - Integral(Integral(w)) + O(z^prec) );
349
z_of_s := S!formal_log(E, prec);
350
sigma_of_s := Evaluate(sigma_of_z, z_of_s);
351
R<s2> := PolynomialRing(RationalField());
352
S<t> := PowerSeriesRing(R);
353
return S!sigma_of_s;
354
end intrinsic;
355
356
intrinsic solve_for_s2(sigma::RngSerElt, p::RngIntElt) -> .
357
{
358
Given a prime p and the formal power series for sigma in terms of s_2,
359
as output by formal_sigma_in_s2, find the element s_2 of Z_p to as high
360
a precision as possible by using that the coefficients of sigma must be
361
362
363
EXAMPLE:
364
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
365
> sig := formal_sigma_in_s2(E, 40);
366
> solve_for_s2(sig, 5);
367
53 + O(5^3)
368
369
3*5^0 + 2*5^2 + O(5^3)
370
}
371
a := 0; // a is what we know about s2 so far.
372
n := 0;
373
S<t> := Parent(sigma);
374
R<s2> := BaseRing(S);
375
prec := AbsolutePrecision(sigma);
376
for i in [3..prec-1] do
377
c := Coefficient(sigma,i);
378
d := Coefficients(c);
379
m := Max([0] cat [Valuation(Denominator(x),p) : x in d]);
380
if m gt 0 then
381
f := p^m*c;
382
// Now viewing f as a poly in s2, it must be 0 modulo p^m, so
383
// as to cancel the denominator.
384
R<x> := PolynomialRing(GF(p));
385
g := R!f;
386
X := Roots(g);
387
assert #X le 1;
388
if #X eq 1 then
389
b := Integers()!X[1][1];
390
a +:= b*p^n;
391
n +:= 1;
392
// Now s2 = a + O(p^(n+1)).
393
z := b+ p*s2;
394
sigma := &+[Evaluate(Coefficient(sigma,r),z)*t^r : r in [0..prec-1]] + O(t^prec);
395
end if;
396
end if;
397
end for;
398
399
end intrinsic;
400
401
402
intrinsic sigma_using_integrality(E::CrvEll, p::RngIntElt,
403
prec::RngIntElt) -> .
404
{
405
Compute the function sigma for E at p using the classical integrality algorithm.
406
407
EXAMPLE:
408
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
409
> sigma_using_integrality(E, 5, 40);
410
O(5^3) + t + O(5^3)*t^2 + (53 + O(5^3))*t^3 - (62 + O(5^3))*t^4 -
411
(23 + O(5^3))*t^5 + (17 + O(5^3))*t^6 - (6 + O(5^2))*t^7 -
412
(58 + O(5^3))*t^8 - (5 + O(5^2))*t^9 + (11 + O(5^2))*t^10 -
413
(2 + O(5))*t^11 + (11 + O(5^2))*t^12 - (1 + O(5))*t^13 +
414
O(5)*t^14 - (1 + O(5))*t^15 + (1 + O(5))*t^16 + O(1)*t^17 +
415
(2 + O(5))*t^18 + O(1)*t^19 + O(1)*t^20 + O(5^-1)*t^21 +
416
O(1)*t^22 + O(5^-1)*t^23 + O(5^-1)*t^24 + O(5^-1)*t^25 +
417
O(5^-1)*t^26 + O(5^-2)*t^27 + O(5^-1)*t^28 + O(5^-2)*t^29 +
418
O(5^-2)*t^30 + O(5^-3)*t^31 + O(5^-2)*t^32 + O(5^-3)*t^33 +
419
O(5^-3)*t^34 + O(5^-3)*t^35 + O(5^-3)*t^36 + O(5^-4)*t^37 +
420
O(5^-3)*t^38 + O(5^-4)*t^39 + O(t^40)
421
}
422
423
sigma := formal_sigma_in_s2(E, prec);
424
s2 := solve_for_s2(sigma, p);
425
R<t> := PowerSeriesRing(Parent(s2));
426
return &+[Evaluate(Coefficient(sigma,n),s2) * t^n : n in [0..prec-1]] + O(t^prec);
427
end intrinsic;
428
429
/**************************************************************************
430
431
Computation of E2(E,omega) and the matrix of absolute Frobenius.
432
433
**************************************************************************/
434
435
intrinsic E2(E::CrvEll, p::RngIntElt, prec::RngIntElt) -> .
436
{
437
Returns value of E2 on the elliptic curve and the matrix
438
of frob on omega=dx/y, eta = x*dx/y, where the elliptic
439
curve is first put in WeiestrassModel form y^2=x^3+ax+b.
440
441
INPUT:
442
E -- elliptic curve
443
p -- prime
444
prec -- p-adic precision for computations.
445
}
446
/*require p ge 5 : "Argument 2 must be at least 5.";*/
447
if p le 3 then
448
print "WARNING: p <= 3.";
449
end if;
450
require IsPrime(p) : "Argument 2 must be prime.";
451
require Conductor(E) mod p ne 0 : "Argument 1 must have good reduction at argument 2.";
452
453
c4, c6 := Explode(cInvariants(MinimalModel(E)));
454
a4 := - c4/(2^4 * 3);
455
a6 := - c6/(2^5 * 3^3);
456
assert IsIsomorphic(E, EllipticCurve([a4,a6]));
457
458
n := 1;
459
d := 2;
460
g := 1;
461
462
/* Set cube = false, so that the following basis is used
463
in the computation of Frob on cohomology:
464
465
omega = dx/y and eta = x(dx/y)
466
467
This is the basis that Katz uses in his computation of E2.
468
*/
469
cube := false;
470
471
N1 := Ceiling((n*g/2)+Log(p,2*Binomial(2*g,g))) + prec;
472
N := N1 + Floor(Log(p,2*N1))+1;
473
474
R1 := quo<Integers(K)|p^N>; // R1 = Z / p^N Z
475
476
L<T> := LaurentSeriesRing(R1);
477
P1 := PolynomialRing(L);
478
X := P1.1;
479
S := quo<P1 | X^3+a4*X+a6 - T^-1>;
480
x := S.1;
481
482
W<X> := PolynomialRing(K);
483
precs := [X^3+a4*X+a6];
484
Append(~precs,Derivative(precs[1]));
485
486
A,B := myXGCD(precs[1],precs[2]);
487
// A,B satisfy A*Q+B*Q'=1 where Q is the lift of poly to char 0.
488
Append(~precs,A);
489
Append(~precs,B);
490
491
/*
492
Compute
493
p*x^(p-1)*(Frob(y))^-1 (or 3 if cube, which isn't true here!)
494
*/
495
496
// print "Computing (y^Frobenius)^(-1)";
497
tyme :=Cputime();
498
x_pow := x^(p-1);
499
poly := PolynomialRing(R1)![a6,a4,0,1];
500
difl := FrobYInv(poly, p, N, x_pow, S,cube) * x_pow;
501
x_pow := x_pow*x;
502
// printf "Expansion time: %o\n",Cputime(tyme);
503
504
/* Calculate the rows of the Frobenius matrix */
505
// print "Reducing differentials modulo cohomology relations.";
506
507
M := MatrixAlgebra(R1,d)!0;
508
i := 1;
509
boundu := p*N + (p div 2) - 1;
510
S1 := PolynomialRing(BaseRing(S));
511
while true do
512
boundl := (p div 2) - Floor((i*p-1)/(d+1));
513
polys,bot := Convert(S1!difl, boundu, boundl, K);
514
diffr := ReduceA([polys[k] : k in [1..Min(1-bot,#polys)]], precs, -bot)+
515
ReduceB([polys[k] : k in [Max(2-bot,1)..#polys]], precs, Max(bot,1), cube);
516
M[i] := RSpace(R1,d)![R1!Coefficient(diffr,k) : k in [0..(d-1)]];
517
if i eq d then
518
break;
519
end if;
520
i +:= 1;
521
difl *:= x_pow;
522
end while;
523
524
/*
525
We know Frob=M to precision N1. Compute the N1-th power of
526
Frob. The second *row* (since MAGMA matrices act from left)
527
contains linear combination of omega and eta that equals
528
Frob^N1(eta).
529
*/
530
A := M^N1;
531
return -12 * (A[2,1]/A[2,2]), M;
532
533
end intrinsic;
534
535
536
/*****************************************************************
537
538
Computing sigma: Cohomology algorithm
539
540
****************************************************************/
541
542
intrinsic sigma_using_e2(E::CrvEll, p::RngIntElt, prec::RngIntElt :
543
e2prec := 50) -> .
544
{
545
Compute and return the formal series for sigma(t) using
546
the Katz/Kedlaya algorithm for computing E2.
547
548
EXAMPLE:
549
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
550
> sigma_using_e2(E, 5, 10);
551
t + O(5^52)*t^2 + (331883312126563673413235306321062928 +
552
O(5^52))*t^3 - (1110223024625156540423631668090820312 +
553
O(5^52))*t^4 - (114496958818289158695517187452183148 +
554
O(5^51))*t^5 + (445053157775380010065972176906892 +
555
O(5^49))*t^6 - (246425046817409275170929836993681 +
556
O(5^48))*t^7 - (248945902282571925665450930262558 +
557
O(5^47))*t^8 - (9730291437202192499870687559126*5 +
558
O(5^47))*t^9 + O(t^10)
559
}
560
prec +:= 2;
561
562
// 1. Compute value of p-adic E_2, using Kedlaya's MW algorithm
563
e2 := E2(E,p,e2prec);
564
565
// 2. Define some notation.
566
567
a1, a2, a3, a4, a6 := Explode(aInvariants(E));
568
S<t> := LaurentSeriesRing(K);
569
570
// 3. Compute formal Weierstrass p function in terms of t.
571
x := S!formal_x(E, prec);
572
wp := x + (a1^2 + 4*a2)/12;
573
574
// 4. Compute wp(z), where z = FormalLog(t)
575
z_of_t := S!formal_log(E, prec);
576
t_of_z := Reverse(z_of_t);
577
wp_of_z := Evaluate(wp, t_of_z);
578
R<z> := S;
579
580
// 5. Compute
581
pr := AbsolutePrecision(wp_of_z);
582
583
// Let g be 1/z^2 - \wp + E2/12. Notice that we compute this
584
// by just ignoring the coefficients of wp_of_z for n=-2.
585
g := e2/12 + &+[z^n*Coefficient(-wp_of_z,n) : n in [-1..pr-1]] + O(z^pr);
586
/* It's *really* weird that it isn't -E2/12, since the canonical
587
eks function is \wp + E2/12.
588
589
There's a factor of -1 in my definition of E2 in the other file,
590
which I think Tate and I got out of Katz's paper. Maybe that
591
sign is wrong? There's some -1 normalization that is different
592
in two papers. */
593
594
// 6. Integrate twice, exponentiate, etc.
595
sigma_of_z := z*my_exp(my_integral(my_integral(g)));
596
sigma_of_t := Evaluate(sigma_of_z, z_of_t);
597
R<t> := PowerSeriesRing(K);
598
sigma := R!sigma_of_t + O(t^(prec-2));
599
return sigma;
600
end intrinsic;
601
602
603
/*****************************************************************
604
605
606
607
****************************************************************/
608
609
intrinsic prep_multiple(P::PtEll, p::RngIntElt) -> PtEll, RngIntElt
610
{
611
Returns an integer n such that n*P is in the subgroup of E that
612
reduces to the identity mod p and lands in the connected component of
613
the Neron model modulo all primes.
614
615
EXAMPLE:
616
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
617
> P := E![0,0];
618
> prep_multiple(P,5);
619
8
620
> prep_multiple(P,11);
621
17
622
}
623
require IsPrime(p) : "Argument 2 must be prime.";
624
625
function prime_divisors(n)
626
return [F[1] : F in Factorization(n)];
627
end function;
628
E := Curve(P);
629
c := LCM([TamagawaNumber(E,q) : q in prime_divisors(Conductor(E))]);
630
N := #ChangeRing(E,GF(p));
631
n := LCM(c,N);
632
return n;
633
end intrinsic;
634
635
636
{
637
Compute and return the value of the unique extension of log to a
638
group homomorphism on all Q_p.
639
640
EXAMPLES:
641
642
> log(K!(1+5+2*5^2));
643
771685135946*5 + O(5^20)
644
> log(K!(5+2*5^2));
645
-159535800608*5 + O(5^20)
646
> log(K!(1/5+2*5^2));
647
-69652953373*5^3 + O(5^20)
648
}
649
p := Prime(Parent(x));
650
u := x * p^(-Valuation(x));
651
return Log(u^(p-1))/(p-1);
652
end intrinsic;
653
654
intrinsic x_coord_d(Q::PtEll) -> RngIntElt
655
{
656
Computes and returns the positive squareroot of the denominator
657
of the x coordinate of Q. Bombs if the denominator is not a perfect square.
658
}
659
t, d := IsSquare(Denominator(Q[1]));
660
assert t;
661
return d;
662
end intrinsic;
663
664
intrinsic my_exp(f::.) -> RngSerLaurElt
665
{
666
Computes exp(f) to the precision of f. This is here because the Exp function
667
built into MAGMA does not work as it should.
668
}
669
return &+[f^n/Factorial(n) : n in [0..AbsolutePrecision(f)]];
670
end intrinsic;
671
672
intrinsic my_is_zero(a::.) -> .
673
{
674
Determines whether or not the p-adic number a is 0.
675
This is here because testing for equality with 0
676
for p-adic numbers in MAGMA does not work as it should.
677
}
678
return Valuation(a) ge Precision(Parent(a));
679
end intrinsic;
680
681
intrinsic my_integral(f::RngSerLaurElt) -> RngSerLaurElt
682
{
683
Integrate Laurent series, without stupid check on coefficient
684
of t^(-1) being 0, which doesn't work in MAGMA, since nonzero
685
for p-adic elements is messed up.
686
}
687
require my_is_zero(Coefficient(f,-1)) : "Coefficient of 1/t must be zero.";
688
pr := AbsolutePrecision(f);
689
t := Parent(f).1;
690
return &+[(Coefficient(f,n)/(n+1))*t^(n+1) : n in [0..pr-1] |
691
not my_is_zero(Coefficient(f,n)) ] + O(t^pr);
692
end intrinsic;
693
694
intrinsic height_function(E::CrvEll, p::RngIntElt, prec::RngIntElt :
695
use_integrality := false) -> .
696
{
697
Returns the CANONICAL p-adic global height function.
698
Note that we normalize this height so it takes values in Z_p, unless
699
p is anomalous in which case it takes values in (1/p)*Z_p.
700
701
INPUT:
702
E -- an elliptic curve over Q,
703
p -- a prime such that E has good ordinary reduction at p,
704
prec -- precision parameter (see below).
705
706
OPTIONS:
707
use_integrality -- this defaults to false.
708
709
* false: Use the E2 algorithm for computing sigma.
710
The heights are computed to precision very
711
close to O(p^prec), but maybe slightly less.
712
I haven't figured out why the precision is
713
slightly less then O(p^prec), since I think
714
it should be exactly O(p^prec), but see that
715
is is by computing the height with larger
716
and larger precision and consider the valuations
717
of the differences.
718
719
* true: Use the integrality of sigma algorithm to
720
compute sigma. The prec argument then determines
721
the number of terms of sigma used in computing s_2.
722
For example, maybe around 40 or 50 terms gives s_2 to
723
precision O(p^3), at least when p=5 and E is 37A.
724
725
EXAMPLES:
726
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
727
> P := E![0,0];
728
> h := height_function(E, 5, 10);
729
> h(P);
730
10120297 + O(5^11)
731
> h := height_function(E, 5, 20);
732
> h(P);
733
-201147061754703 + O(5^21)
734
> Valuation(-201147061754703 - 10120297, 5);
735
8
736
> h := height_function(E, 5, 30);
737
> h(P);
738
-1305176752965909410953 + O(5^31)
739
> Valuation(-1305176752965909410953 + 201147061754703, 5);
740
18
741
> h := height_function(E, 5, 60 : use_integrality := true); // slow
742
> h(P);
743
> -3 + O(5^2) // pitiful precision
744
}
745
require Conductor(E) mod p ne 0 and ap(E,p) mod p ne 0 : "Curve must have good ordinary reduction at p.";
746
747
if use_integrality then
748
sig := sigma_using_integrality(E, p, prec) ;
749
else
750
sig := sigma_using_e2(E, p, prec : e2prec := prec) ;
751
end if;
752
function h(P)
753
assert Curve(P) eq E;
754
n := prep_multiple(P, p);
755
Q := n*P;
756
d := x_coord_d(Q); // positive sqrt of denominator
757
t := t_val(Q);
758
v := Evaluate(sig,t);
759
HQ := log(v/d);
760
HP := HQ/n^2;
761
return HP/p; // note the normalization
762
end function;
763
return h;
764
end intrinsic;
765
766
intrinsic regulator(E::CrvEll, p::RngIntElt, prec::RngIntElt) -> .
767
{
768
Computes and returns the p-adic regulator of E, which is the
769
discriminant of the height pairing on the Mordell-Weil group E(Q).
770
771
INPUT:
772
E -- an elliptic curve over Q.
773
p -- a good ordinary prime for E.
774
prec -- the precision that is input to the height_function command, which
775
is computed using the E2 algorithm.
776
777
EXAMPLES:
778
> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
779
> regulator(E,5,20);
780
138360994885922 + O(5^21)
781
> regulator(E,7,20);
782
4709403600911866 + O(7^20)
783
784
> regulator(E,7,20);
785
22975764581280320 + O(7^20)
786
}
787
require Conductor(E) mod p ne 0 and ap(E,p) mod p ne 0 :
788
"Curve must have good ordinary reduction at p.";
789
790
h := height_function(E, p, prec);
791
G, f := MordellWeilGroup(E);
792
B := [f(G.i) : i in [1..Ngens(G)] | Order(G.i) eq 0];
793
if #B eq 0 then
794
return 1;
795
end if;
796
function pairing(P,Q)
797
return (h(P+Q)-h(P)-h(Q))/2;
798
end function;
799
K := Parent(h(B[1]));
800
M := MatrixAlgebra(K,#B)!0;
801
for i in [1..#B] do
802
for j in [i..#B] do
803
aij := pairing(B[i],B[j]);
804
M[i,j] := aij;
805
M[j,i] := aij;
806
end for;
807
end for;
808
return Determinant(M);
809
end intrinsic;
810
811
812
/***********************************************************************
813
814
Miscellaneous functions that are useful to have around.
815
816
***********************************************************************/
817
intrinsic ap(E::CrvEll, p::RngIntElt) ->RngIntElt
818
{
819
Returns #E(F_p) - (p + 1).
820
}
821
return TraceOfFrobenius(ChangeRing(MinimalModel(E),GF(p)));
822
end intrinsic;
823
824
intrinsic good_ordinary_primes(E::CrvEll, pmax::RngIntElt) -> SeqEnum
825
{
826
Compute and return the primes p>=5 that are good ordinary for E and are
827
less than pmax.
828
}
829
N := Conductor(E);
830
return [p : p in [5..pmax] | IsPrime(p) and
831
N mod p ne 0 and
832
ap(E,p) mod p ne 0];
833
end intrinsic;
834
835
836
{
837
I can't figure out how to print p-adics correctly in MAGMA (as power
838
series in p) anymore. This function isn't very refined though, since
839
e.g., it outputs 3*5^0 + 2*5^2 instead of 3+2*5^2.
840
841
EXAMPLE:
842
843
> a := K!9484945; a;
844
1896989*5 + O(5^21)
845
846
4*5^1 + 2*5^2 + 4*5^3 + 2*5^6 + 5^7 + 4*5^8 + 4*5^9 + O(5^10)
847
}
848
z := Integers()!x;
849
p := Prime(Parent(x));
850
i := 0;
851
s := "";
852
while z ne 0 and i lt Precision(Parent(x)) do
853
c := z mod p;
854
if c ne 0 then
855
if c ne 1 then
856
s *:= Sprintf("%o*", c);
857
end if;
858
s *:= Sprintf("%o^%o + ", p, i);
859
end if;
860
z := z div p;
861
i := i + 1;
862
end while;
863
s *:= Sprintf("O(%o^%o)", p, i);
864
return s;
865
end intrinsic;
866
867
868
/* A computation for Mazur-Rubin-Stein:
869
870
Email from Mazur:
871
872
of X_0(11) of Mordell-Weil rank 2 (p=3)
873
From: barry mazur <mazur@math.harvard.edu>
874
To: "William A. Stein" <was@math.harvard.edu>
875
CC: Karl Rubin <krubin@math.uci.edu>
876
Date: 2004-11-29 09:58 am
877
878
Hi William,
879
880
Consider E = X_0(11) and look for quadratic characters \chi of
881
conductor d=d(\chi) such that \chi(-11) = +1, and such that the L
882
function L(E,\chi,s) vanishes at the central point s=1. Then we can
883
generally expect that the Mordell-Weil rank of E_{\chi} is 2
884
(although, of course, very very occasionally it might be \ge 4). The
885
small values of d=d(\chi) where \chi(-11) = +1, and L(E,\chi,1) = 0
886
are, I think,
887
d = -47, -103, -119, -159, -344, -455, -488 (for negative
888
conductors)
889
d = 232, 265, 273, 344, 364, 401, 421, 476, 488 (for positive
890
conductors)
891
892
We want to fix a prime p of good ordinary reduction for E_{\chi}, and
893
compute to some accuracy the 2x2 matrix giving the p-adic height
894
pairing relative to some basis of the Mordell-Weil group of E_{\chi}.
895
In my article with Karl, we show that the Selmer \Lambda-module of
896
E_{\chi} is pretty much understood if we know this matrix. Actually,
897
we aren't interested in the result UNLESS this matrix has determinant
898
divisible by p. So, keep to p=3, (or perhaps also p=7,13,17,23), and
899
consider only the d's that are prime to p. Of course, one might also
900
do a triage by:
901
902
1) We might just try to compute the second derivative mod p of the
903
p-adic L function L_p(E,\chi,T) where T = (1+p)^s as usual, in hopes
904
of getting zero mod p. Only if so, we could:
905
906
2) Look for two generators of Mordell-Weil, and try to compute the
907
908
909
What do you think?
910
*/
911
912
intrinsic mazur_comp(d::RngIntElt)
913
{}
914
915
916
time G,f := MordellWeilGroup(F);
917
if Rank(F) eq 0 then
918
print "MAGMA couldn't find points for d=",d, " so skipping.";
919
return;
920
end if;
921
// Compute the good ordinary primes up to 100.
922
P := good_ordinary_primes(F, 30);
923
// For each, compute the regulator to 20 terms of precision.
924
for p in P do
925
reg := regulator(F, p, 20);
926
print <d, p, reg, Valuation(reg)>;
927
end for;
928
end intrinsic;
929