Author: William A. Stein
Compute Environment: Ubuntu 18.04 (Deprecated)
1/*****************************************************************
2
4                  cyclotomic height pairing.
5
6 (c) 2004, William Stein, was@math.harvard.edu
7
8 ****************************************************************/
9
10/*****************************************************************
111. 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
24import "kedlaya.m": myXGCD, FrobYInv, Convert, ReduceA, ReduceB;
25
26/*****************************************************************
27   2. How do you use this package?
28
29First try the following calculations and see if they work:
30
31> Attach("kedlaya.m");
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
39Also, compute some regulators:
40
41> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
42> regulator(E,5,20);
43138360994885922 + O(5^21)
44> regulator(E,7,20);
454709403600911866 + O(7^20)
47> regulator(E,7,20);
4822975764581280320 + O(7^20)
49
50Now you should know enough to be able to make use of the package.
51However, you should probably read through the source code below, To
52make this easier, I've put a description of *every* function with a
53complete example of usage.  I've also grouped the functions into
54logical sections.
55
56****************************************************************/
57
58
59/*****************************************************************
60
61  Formal groups of elliptic curves.
62
63 ****************************************************************/
64
65R<t> := LaurentSeriesRing(RationalField());
66
67intrinsic formal_w(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
68{
69Compute and return the formal power series w(t) defined in Chapter IV, Section 1
70of Silverman's AEC book.  This function is useful because Silverman gives
71expressions in w(t) for other formal power series related to formal groups.
72
73EXAMPLE:
74> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
75> formal_w(E,19);
76t^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;
90end intrinsic;
91
92intrinsic formal_x(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
93{
94Compute and return the formal group power series defined by the
95function x on E.
96
97EXAMPLES:
98> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
99> formal_x(E,15);
100t^-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);
104end intrinsic;
105
106intrinsic formal_y(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
107{
108Compute and return the formal group power series defined by the
109function y on E.
110
111EXAMPLES:
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);
118end intrinsic;
119
120
121intrinsic formal_omega(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
122{
123Compute and return the formal differential in terms of t, i.e., return
124the function f(t) such that f(t)dt is the formal differential omega.
125
126ALGORITHM:
127    Use that omega = dx/(2*y + a1*x + a3) = f(t) dt.
128
129EXAMPLES:
130> formal_omega(E,10);
1311 + 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
141end intrinsic;
142
143intrinsic formal_log(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
144{
145The formal log in terms of t.  This is sometimes written z=z(s).  It is the
146integral with respect to t of the formal differential.
147
148EXAMPLE:
149> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
150> formal_log(E,10);
151t + 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));
158end intrinsic;
159
160intrinsic formal_inverse(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
161{
162Compute the power series in Z[[t]] corresponding to the inverse
163map in the formal group on E.
164
165EXAMPLES:
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);
174end intrinsic;
175
176intrinsic formal_group_law(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
177{
178Compute and return the formal group law on E as an element of Z[[t,t2]].
179The precision must be at least 4.
180
181ALGORITHM:
182  Use the formula at the end of Section 1 of Chapter 4 of Silverman AEC,
183on page 114.  This is the formula for "z_3" there.  Note though that
184the FORMULA IS *WRONG* in Silverman's book!!  The a_1 and a_3 in that
185formula should have minus signs in front, but they don't.
186
187
188EXAMPLE:
189> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
190> formal_group_law(E,4);
191t + 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);
211end intrinsic;
212
213intrinsic formal_mult(E::CrvEll, n::RngIntElt, prec::RngIntElt) -> RngSerLaurElt
214{
215Compute and return the formal power series corresponding to multiplication by n on E,
216to precision prec.
217
218EXAMPLE:
219> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
220>  formal_mult(E,2,5);
2212*t - 7*t^4 + 12*t^5 + 4*t^7 + O(t^8)
222>  formal_mult(E,3,5);
2233*t - 39*t^4 + 96*t^5 + 234*t^7 + O(t^8)
224>  formal_mult(E,7,5);
2257*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);
242end intrinsic;
243
244intrinsic formal_divpoly(E::CrvEll, m::RngIntElt, prec::RngIntElt) -> RngSerElt
245{
246Compute and return the formal division "polynomial" f_m for E to
247precision prec.
248
249EXAMPLES:
250> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
251> formal_divpoly(E,2,5);
2522*t^-3 - 3 + 2*t - 2*t^3 + O(t^4)
253> formal_divpoly(E,3,5);
2543*t^-8 - 12*t^-5 + 6*t^-4 + 9*t^-2 + O(t^-1)
255> formal_divpoly(E,7,5);
2567*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));
269end intrinsic;
270
271intrinsic weierstrass_p(E::CrvEll, prec::RngIntElt) -> FldPrElt
272{
273   Compute and return the Weierstrass p-series for E to precision prec.
274
275EXAMPLES:
276> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
277> weierstrass_p(E,9);
278t^-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)));
294end intrinsic;
295
296
297intrinsic t_val(Q::PtEll)  -> FldRatElt
298{
299Given a point (x,y) on an elliptic curve, return t = -x/y.  (This function doesn't
300do much!)
301
302EXAMPLE:
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);
3122/5
313}
314   require Q ne 0 : "The denominator y must be nonzero.";
315   return -Q/Q;
316end intrinsic;
317
318
319/*****************************************************************
320
321    Computing sigma: Integrality algorithm
322
323 ****************************************************************/
324
325intrinsic formal_sigma_in_s2(E::CrvEll, prec::RngIntElt) -> RngSerElt
326{
327Compute a formal power series for sigma with coefficient polynomials
328in the coefficient s_2 of t^3.
329
330NOTE: This is adapted from Christian Wuthrich's PARI code, which is also the
331algorithm Mazur-Tate and I found on 2004-06-09.
332
333EXAMPLE:
334> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
335> formal_sigma_in_s2(E, 10);
336t + 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;
354end intrinsic;
355
356intrinsic solve_for_s2(sigma::RngSerElt, p::RngIntElt) -> .
357{
358Given a prime p and the formal power series for sigma in terms of s_2,
359as output by formal_sigma_in_s2, find the element s_2 of Z_p to as high
360a precision as possible by using that the coefficients of sigma must be
362
363EXAMPLE:
364> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
365> sig := formal_sigma_in_s2(E, 40);
366> solve_for_s2(sig, 5);
36753 + O(5^3)
3693*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( 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;
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;
399end intrinsic;
400
401
402intrinsic sigma_using_integrality(E::CrvEll, p::RngIntElt,
403				  prec::RngIntElt) -> .
404{
405Compute the function sigma for E at p using the classical integrality algorithm.
406
407EXAMPLE:
408> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
409> sigma_using_integrality(E, 5, 40);
410O(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);
427end intrinsic;
428
429/**************************************************************************
430
431 Computation of E2(E,omega) and the matrix of absolute Frobenius.
432
433 **************************************************************************/
434
435intrinsic E2(E::CrvEll, p::RngIntElt, prec::RngIntElt) -> .
436{
437Returns value of E2 on the elliptic curve and the matrix
438of frob on omega=dx/y, eta = x*dx/y, where the elliptic
439curve is first put in WeiestrassModel form y^2=x^3+ax+b.
440
441INPUT:
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;
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));
485
486   A,B := myXGCD(precs,precs);
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.";
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
533end intrinsic;
534
535
536/*****************************************************************
537
538    Computing sigma: Cohomology algorithm
539
540 ****************************************************************/
541
542intrinsic 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
548EXAMPLE:
549> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
550> sigma_using_e2(E, 5, 10);
551t + 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.
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;
600end intrinsic;
601
602
603/*****************************************************************
604
605    Computing the p-adic Height Function
606
607 ****************************************************************/
608
609intrinsic prep_multiple(P::PtEll, p::RngIntElt) -> PtEll, RngIntElt
610{
611Returns an integer n such that n*P is in the subgroup of E that
612reduces to the identity mod p and lands in the connected component of
613the Neron model modulo all primes.
614
615EXAMPLE:
616> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
617> P := E![0,0];
618> prep_multiple(P,5);
6198
620> prep_multiple(P,11);
62117
622}
623   require IsPrime(p) : "Argument 2 must be prime.";
624
625   function prime_divisors(n)
626      return [F : 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;
633end intrinsic;
634
636{
637Compute and return the value of the unique extension of log to a
638group homomorphism on all Q_p.
639
640EXAMPLES:
642> log(K!(1+5+2*5^2));
643771685135946*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);
652end intrinsic;
653
654intrinsic x_coord_d(Q::PtEll) -> RngIntElt
655{
656Computes and returns the positive squareroot of the denominator
657of the x coordinate of Q.  Bombs if the denominator is not a perfect square.
658}
659   t, d := IsSquare(Denominator(Q));
660   assert t;
661   return d;
662end intrinsic;
663
664intrinsic my_exp(f::.) -> RngSerLaurElt
665{
666Computes exp(f) to the precision of f.  This is here because the Exp function
667built into MAGMA does not work as it should.
668}
669   return &+[f^n/Factorial(n) : n in [0..AbsolutePrecision(f)]];
670end intrinsic;
671
672intrinsic my_is_zero(a::.) -> .
673{
674Determines whether or not the p-adic number a is 0.
675This is here because testing for equality with 0
676for p-adic numbers in MAGMA does not work as it should.
677}
678   return Valuation(a) ge Precision(Parent(a));
679end intrinsic;
680
681intrinsic 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);
692end intrinsic;
693
694intrinsic height_function(E::CrvEll, p::RngIntElt, prec::RngIntElt :
695            use_integrality := false) -> .
696{
697Returns the CANONICAL p-adic global height function.
698Note that we normalize this height so it takes values in Z_p, unless
699p is anomalous in which case it takes values in (1/p)*Z_p.
700
701INPUT:
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
706OPTIONS:
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
725EXAMPLES:
726> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
727> P := E![0,0];
728> h := height_function(E, 5, 10);
729> h(P);
73010120297 + O(5^11)
731> h := height_function(E, 5, 20);
732> h(P);
733-201147061754703 + O(5^21)
734> Valuation(-201147061754703 - 10120297, 5);
7358
736> h := height_function(E, 5, 30);
737> h(P);
738-1305176752965909410953 + O(5^31)
739> Valuation(-1305176752965909410953  + 201147061754703, 5);
74018
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;
764end intrinsic;
765
766intrinsic 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
771INPUT:
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
777EXAMPLES:
778> E := EllipticCurve([ 0, 0, 1, -1, 0 ]);
779> regulator(E,5,20);
780138360994885922 + O(5^21)
781> regulator(E,7,20);
7824709403600911866 + O(7^20)
784> regulator(E,7,20);
78522975764581280320 + 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));
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);
809end intrinsic;
810
811
812/***********************************************************************
813
814   Miscellaneous functions that are useful to have around.
815
816 ***********************************************************************/
817intrinsic ap(E::CrvEll, p::RngIntElt) ->RngIntElt
818{
819  Returns #E(F_p) - (p + 1).
820}
821   return TraceOfFrobenius(ChangeRing(MinimalModel(E),GF(p)));
822end intrinsic;
823
824intrinsic 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];
833end intrinsic;
834
836{
837I can't figure out how to print p-adics correctly in MAGMA (as power
838series in p) anymore.   This function isn't very refined though, since
839e.g., it outputs 3*5^0 + 2*5^2 instead of 3+2*5^2.
840
841EXAMPLE:
843> a := K!9484945; a;
8441896989*5 + O(5^21)
8464*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;
865end intrinsic;
866
867
868/* A computation for Mazur-Rubin-Stein:
869
870Email from Mazur:
872    of X_0(11) of Mordell-Weil rank 2  (p=3)
873From: barry mazur <mazur@math.harvard.edu>
874To: "William A. Stein" <was@math.harvard.edu>
875CC: Karl Rubin <krubin@math.uci.edu>
876Date: 2004-11-29 09:58 am
877
878  Hi William,
879
880Consider  E = X_0(11) and look for quadratic characters \chi of
881conductor d=d(\chi)  such that \chi(-11) = +1, and such that the L
882function L(E,\chi,s) vanishes at the central point s=1.  Then we can
883generally expect that the Mordell-Weil rank of E_{\chi} is 2
884(although, of course, very very occasionally it might be \ge 4).   The
885small values of d=d(\chi) where  \chi(-11) = +1, and  L(E,\chi,1) = 0
886are, I think,
887   d = -47, -103, -119, -159, -344, -455, -488      (for negative
888conductors)
889   d = 232, 265, 273, 344, 364, 401, 421, 476, 488  (for positive
890conductors)
891
892We want to fix a prime p of good ordinary reduction for  E_{\chi}, and
893compute to some accuracy the 2x2 matrix giving the p-adic height
894pairing relative to some basis of the Mordell-Weil group of E_{\chi}.
895In my article with Karl, we show that the Selmer \Lambda-module of
896E_{\chi} is pretty much understood if we know this matrix.  Actually,
897we aren't interested in the result UNLESS this matrix has determinant
898divisible by p.  So, keep to  p=3, (or perhaps also p=7,13,17,23), and
899consider only the d's that are prime to p.  Of course, one might also
900do a triage by:
901
9021) We might just try to compute the second derivative mod p of the
903p-adic L function L_p(E,\chi,T)  where T = (1+p)^s as usual, in hopes
904of getting zero mod p. Only if so, we could:
905
9062) Look for two generators of Mordell-Weil, and try to compute the
908
909What do you think?
910*/
911
912intrinsic mazur_comp(d::RngIntElt)
913{}
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;
928end intrinsic;
929