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