CoCalc Public Fileswww / talks / harvard-talk-2004-12-08 / e2heights / formal_orig.m
Author: William A. Stein
Compute Environment: Ubuntu 18.04 (Deprecated)
1R<t> := LaurentSeriesRing(RationalField());
2Q := RationalField();
3
4intrinsic WeierstrassP(E::CrvEll, prec::RngIntElt) -> FldPrElt
5{}
6   vprint Heegner :  "Computing WeierstrassP";
7   L<w> := PolynomialRing(Rationals());
8   R<t> := LaurentSeriesRing(L);
9   p := t^(-2);
10   c4, c6 := Explode(cInvariants(E));
11   a := 4; b := -c4/12; c := -c6/216;
12   for n in [1..prec] do
13      ptest := p + w*t^n + O(t^(n+1));
14      ptest_prime := Derivative(ptest);
15      f := ptest_prime^2 - a*ptest^3 - b*ptest - c;
16      x := Coefficient(f,-4+n);
17      coeff := -Coefficient(x,0)/Coefficient(x,1);
18      p := p + coeff*t^n;
19   end for;
20   vprint Heegner :  "Done computing Weierstrass P";
21   return LaurentSeriesRing(Rationals())!(p + O(t^(prec+1)));
22end intrinsic;
23
24intrinsic formal_w(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
25{}
26   w := t^3 + O(t^prec);
27   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
28   for n in [4..prec-1] do
29      w := w +  t^n*(a1*Coefficient(w,n-1) +
30                a2*Coefficient(w,n-2) +
31                a3*&+[Integers()|Coefficient(w,i)*Coefficient(w,n-i) : i in [3..n-3]] +
32                a4*&+[Integers()|Coefficient(w,i)*Coefficient(w,n-1-i) : i in [3..n-4]] +
33                a6*&+[Integers()|Coefficient(w,i)*Coefficient(w,j)*Coefficient(w,n-i-j) :
34                         i in [3..n-3], j in [3..n-3] | n-i-j ge 3]);
35   end for;
36   return w;
37end intrinsic;
38
39intrinsic formal_x(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
40{}
41   return t/formal_w(E,prec+5);
42end intrinsic;
43
44intrinsic formal_y(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
45{}
46   return -1/formal_w(E,prec+6);
47end intrinsic;
48
49intrinsic formal_omega(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
50{Compute the formal differential in terms of t, i.e., return the function f(t)
51such that f(t)dt is the formal differerntial omega.
52
53Use that omega = dx/(2*y + a1*x + a3) = f(t) dt.
54}
55
56   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
57   x := formal_x(E,prec);
58   y := formal_y(E,prec);
59   return Derivative(x)/(2*y + a1*x + a3);
60end intrinsic;
61
62intrinsic formal_log(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
63{The formal log in terms of t.  This is sometimes written z=z(s).  It is the
64integral with respect to t of the formal differential.}
65   y := formal_y(E,prec);
66   x := formal_x(E,prec);
67   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
68   return Integral(1/(2*y + a1*x + a3) * Derivative(x));
69end intrinsic;
70
71intrinsic formal_inverse(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
72{}
73   x := formal_x(E,prec);
74   y := formal_y(E,prec);
75   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
76   return x/(y + a1*x + a3);
77end intrinsic;
78
79intrinsic formal_group(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
80{}
81   t1 := t;
82   S<t2> := LaurentSeriesRing(R);
83   w := formal_w(E,prec);
84   function tsum(n)
85      return &+[t2^m * t1^(n-m-1) : m in [0..n-1]];
86   end function;
87   lam := &+[tsum(n)*Coefficient(w,n) : n in [3..prec-1]];
88
89//   lam  := (Evaluate(w,t2+O(t2^prec)) - Evaluate(w,t1+O(t1^prec)))/(t2-t1);
90
91   w1 := Evaluate(w,t1+O(t^prec));
92   nu := w1 - lam*(t1+O(t^prec));
93   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
94/*   t3 := -t1 - t2 + \
95         (a1*lam + a3*lam^2 - a2*nu - 2*a4*lam*nu - 3*a6*lam^2*nu)/
96         (1 + a2*lam + a4*lam^2 + a6*lam^3) + O(t1^prec) + O(t2^prec);
97*/
98   t3 := -t1 - t2 + \
99         (-a1*lam - a3*lam^2 - a2*nu - 2*a4*lam*nu - 3*a6*lam^2*nu)/
100         (1 + a2*lam + a4*lam^2 + a6*lam^3) + O(t1^prec) + O(t2^prec);
101   inv := formal_inverse(E, prec);
102   return Evaluate(inv, t3);
103end intrinsic;
104
105
106
107intrinsic formal_mult(E::CrvEll, n::RngIntElt, prec::RngIntElt) -> RngSerLaurElt
108{}
109   if n eq -1 then
110      return formal_inverse(E,prec);
111   end if;
112
113   // The following is less direct, but it works more quickly...
114   L := FunctionField(E);
115   EL := BaseExtend(E,L);
116   Q := EL![L.1,L.2];
117   x,y := Explode(Eltseq(n*Q));
118   a := -x/y;
119   xx := formal_x(E,prec);
120   yy := formal_y(E,prec);
121   phi := hom<Parent(x) -> Parent(xx) | xx,yy>;
122   return phi(a);
123
124// more direct method... (not used if the above isn't commented out)
125
126   if n eq 1 then
127      return t;
128   end if;
129   if n eq -1 then
130      return formal_inverse(E,prec);
131   end if;
132   if n lt 0 then
133      return Evaluate(formal_inverse(E,prec),formal_mult(E,-n,prec));
134   end if;
135   // Now n is positive and >= 2.
136   F := formal_group(E,prec);
137   g := t;  // 1
138   for m in [2..n] do
139      g := Evaluate(F,g);
140   end for;
141   return g;
142
143/*
144   // Use a binary powering algorithm to compute g = [n](t).
145   g := t + O(t^prec);  // [1]
146   ans := 0;
147   while n gt 0 do
148      if n mod 2 ne 0 then
149         ans := Evaluate(Evaluate(F,ans), g);
150      end if;
151      g := Evaluate(Evaluate(F,g),g);   // [m] <-- [2m], where g(t) <--> [m].
152      n := n div 2;
153   end while;
154   return ans;
155*/
156end intrinsic;
157
158intrinsic indeterminate_sigma(prec::RngIntElt) -> RngSerElt
159{}
160   assert prec ge 3;
161   S := PolynomialRing(Q, prec-2);
162   AssignNames(~S,[Sprintf("s%o",i) : i in [1..prec-2]]);
163   T<t> := LaurentSeriesRing(S);
164   // the 2 ------------------------> is for special case s1=0.
165   //return t + &+[ S.i*t^(i+1) : i in [2..prec-2]] + O(t^prec);
166   return t + &+[ S.i*t^(i+1) : i in [1..prec-2]] + O(t^prec);
167end intrinsic;
168
169intrinsic formal_divpoly(E::CrvEll, m::RngIntElt, prec::RngIntElt) -> RngSerElt
170{}
171   if m eq 1 then
172      return R!1;
173   end if;
174   if m eq -1 then
175      return R!(-1);
176   end if;
177   if m eq 2 then
178      return Sqrt(Evaluate(DivisionPolynomial(E,m),formal_x(E,prec)));
179   end if;
180   assert IsOdd(m);
181   return Evaluate(DivisionPolynomial(E,m),formal_x(E,prec));
182end intrinsic;
183
184intrinsic sigma_relation(E::CrvEll, m::RngIntElt,
185                  numvars::RngIntElt, prec::RngIntElt) -> SeqEnum
186{}
187   sig := indeterminate_sigma(numvars+2);
188   T<t> := Parent(sig);
189   mult := T!Eltseq(formal_mult(E,m,prec)) * t;
190   fm := T!Eltseq(formal_divpoly(E,m,prec)) * t^(1-m^2);
191   return (t^(m^2-1)*Evaluate(sig,mult) - t^(m^2-1)*sig^(m^2)*fm);
192end intrinsic;
193
194/* The two functions formal_sigma_in_s2 and formal_sigma_in_s2_2 compute
195   the same thing in totally different ways.  The first one is more
196   direct and hence vastly faster.
197*/
198
199intrinsic formal_sigma_in_s2(E::CrvEll, prec::RngIntElt) -> RngSerElt
200{}
201   // This is adapted from Christian's PARI code, which is also exactly
202   //  what Mazur-Tate worked out on 2004-06-09.
203
204   R<s2> := FieldOfFractions(PolynomialRing(RationalField()));
205   S<z> := LaurentSeriesRing(R);
206   w := S!WeierstrassP(E,prec);
207   w := &+[Coefficient(w,n)*z^n : n in [0..prec-1]] + O(z^prec);
208   sigma_of_z := z*Exp( s2*z^2 - Integral(Integral(w)) + O(z^prec) );
209   z_of_s := S!formal_log(E, prec);
210   sigma_of_s := Evaluate(sigma_of_z, z_of_s);
211   R<s2> := PolynomialRing(RationalField());
212   S<t> := PowerSeriesRing(R);
213   return S!sigma_of_s;
214end intrinsic;
215
216intrinsic formal_sigma_in_s2_2(E::CrvEll, prec::RngIntElt) -> RngSerElt
217{}
218   m := 2;
219   rels := Eltseq(sigma_relation(E,m,prec,prec+2));
220   S := Parent(rels[1]);
221   i := 1;
222   seen := [2];
223   numvars := prec;
224   known := [S.i : i in [1..numvars]];
225   for r in rels do
226      mon := Monomials(r);
227      cof := Coefficients(r);
228      z := 0;
229      for j in [1..#mon] do
230         for i in [1..numvars] do
231            if mon[j] eq S.i and not (i in seen) then
232               z := &+[S|-mon[k]*cof[k]/cof[j] : k in [1..#mon] | k ne j];
233               //print "\nfor var s",i," have expression\n",z;
234               phi := hom<S -> S | known>;
235               val := phi(z);
236               //printf "\ns[%o] := %o;\n", i, val;
237               //print Valuation(LCM([1] cat [Denominator(x) : x in Coefficients(val)]) , 5);
238               known[i] := val;
239               Append(~seen,i);
240            end if;
241         end for;
242      end for;
243   end for;
244   R<s2> := PolynomialRing(RationalField());
245   z := [R|0 : i in [1..numvars]];
246   z[2] := s2;
247   phi := hom<S -> R | z>;
248   S<t> := PowerSeriesRing(R);
249   return S!([0,1] cat [phi(c) : c in known]) + O(t^(#known+1));
250end intrinsic;
251
252
253
254/*******************************************************
255
256Karl Rubin and I were wondering whether you could do a p-adic height
257computation, the point of which is to figure out whether the p-adic
258regulator of elliptic curves has a "tendency" to be not divisible by p. We
259want to avoid a few degenerate cases, so the smallest p that is useful to
260us is p=5. Is the strategy-- described below--, that systematically runs
262(non-3-Eisenstein) elliptic curve of conductor 37 reasonably easy to do?
263
264
265%%%%%%%%%%%%%%%%%%%
266Our untwisted equation is E:  y^2 = 4x^3-4x+1. For any nonzero integer D we
267consider the twisted elliptic curve  E_D:  D.y^2 = 4x^3-4x+1.  Put F(x): =
2684x^3-4x+1. Let x run through integers
269
270  a=  0, \pm 1, \pm 2, \pm3, ...
271
272and for each of these a's,
273
2741) Factor F(a) = d.b^2,  where d is square-free;
275
2762) throw out the ones where d is divisible by 5, or where the parity of the
277rank of E_d is even;
278
2793) for the rest, viewing P:=(a,b) as a rational point-- it will be of
280infinite order--on the elliptic curve E_d, so one can compute the 5-adic
281height of P on the elliptic curve E_d: this should be a 5-adic integer if
282the height is normalized right, and there should be a good number of
283instances where  this height is a 5-adic unit; throw them out, and
284
2854) print the rest.
286%%%%%%%%%%%%%%%%
287
288Note that we only want to know the 5-adic height modulo 5.
289
290Then it would be our job to figure out whether strange things occur for
291this printed list of instances: e.g., the point P might be divisible by 5,
292or the Mordell-Weil rank might be > 2, or other stuff..., and we want to
293know what other stuff there might be.
294
295****************************************************/
296
297intrinsic point(a::RngIntElt) -> PtEll, RngIntElt, BoolElt
298{Returns:
299     * point P_d on a minimal model for E_d,
300     * the integer d,
301     * a boolean -- true <==> an_rank(E_d) is odd.
302}
303   d, b := SquareFree(4*a^3-4*a+1);
304   E := EllipticCurve([-16*d^2, 16*d^3]);
305   P := E![4*d*a, 4*d^2*b];
306   EM, phi := MinimalModel(E);
307   P := phi(P);
308   return P, d, RootNumber(E) eq 1;
309end intrinsic;
310
311function prime_divisors(n)
312   return [F[1] : F in Factorization(n)];
313end function;
314
315intrinsic prep_multiple(P::PtEll, p::RngIntElt) -> PtEll, RngIntElt
316{An integer n such that n*P is in the subgroup of
317E that reduces to the identity mod p and lands in the connected
318component of the Neron model modulo all primes.}
319   require IsPrime(p) : "Argument 2 must be prime.";
320
321   E := Curve(P);
322   c := LCM([TamagawaNumber(E,q) : q in prime_divisors(Conductor(E))]);
323   N := #ChangeRing(E,GF(p));
324   n := LCM(c,N);
325   return n;
326end intrinsic;
327
329{}
330   p := Prime(Parent(x));
331   u := x * p^(-Valuation(x));
332   return Log(u^(p-1))/(p-1);
333end intrinsic;
334
335intrinsic x_coord_d(Q::PtEll) -> RngIntElt
336{}
337   t, d := IsSquare(Denominator(Q[1]));
338   assert t;
339   return d;
340end intrinsic;
341
342intrinsic t_val(Q::PtEll)  -> FldRatElt
343{Given input (x,y), return t = -x/y.}
344   return -Q[1]/Q[2];
345end intrinsic;
346
347intrinsic formal_val_of_sigma_in_s2(E::CrvEll, t::FldRatElt,
348                     numvars::RngIntElt, prec::RngIntElt) -> .
349{}
350   sigma := formal_sigma_in_s2(E, numvars, prec);
351   return Evaluate(sigma, t);
352end intrinsic;
353
354intrinsic mult_height_of_associated_prepared_point(P::PtEll, p::RngIntElt) -> .
355{Returns multiplicative p-adic height as series in sigma.}
356   n := prep_multiple(P, p);
357   Q := n*P;
358   d := x_coord_d(Q);  // positive sqrt of denominator
359   HQ := formal_val_of_sigma_in_s2(Curve(P), t_val(Q), 5, 10)/d;
360   return HQ, n;
361end intrinsic;
362
363intrinsic solve_for_s2(sigma::RngSerElt, p::RngIntElt) -> .
364{}
365   a := 0;   // a is what we know about s2 so far.
366   n := 0;
367   S<t> := Parent(sigma);
368   R<s2> := BaseRing(S);
369   prec := AbsolutePrecision(sigma);
370   for i in [3..prec-1] do
371      c := Coefficient(sigma,i);
372      d := Coefficients(c);
373      m := Max([0] cat [Valuation(Denominator(x),p) : x in d]);
374      if m gt 0 then
375         f := p^m*c;
376         // Now viewing f as a poly in s2, it must be 0 modulo p^m, so
377         // as to cancel the denominator.
378         R<x> := PolynomialRing(GF(p));
379         g := R!f;
380         X := Roots(g);
381         assert #X le 1;
382         if #X eq 1 then
383            b := Integers()!X[1][1];
384            a +:= b*p^n;
385            n +:= 1;
386            // Now s2 = a + O(p^(n+1)).
387            z := b+ p*s2;
388            sigma := &+[Evaluate(Coefficient(sigma,r),z)*t^r : r in [0..prec-1]] + O(t^prec);
389         end if;
390      end if;
391   end for;
393end intrinsic;
394
395intrinsic compute_s2(E::CrvEll, p::RngIntElt, prec::RngIntElt, s1::.) -> RngElt
396{
397  Compute the quantity s2, which is the coefficient of t^3 in the expansion
398  of the p-adic sigma function associated to E.
399}
400/*
401  If I computed everything correctly, we should have that
402
403     s2 = -(a1^2+4*a2+E2(E))/24 - (a2 - a1*s1 - s1^2)/2
404
405*/
406    e2 := E2(E,p,prec);
407    a1, a2, a3, a4, a6 := Explode(aInvariants(E));
408    return (-(e2 + a1^2 + 4*a2)/12 -a2 +a1*s1 + s1^2)/2;
409end intrinsic;
410
411intrinsic val_of_sigma(Q::PtEll, p::RngIntElt,
412                       prec::RngIntElt) -> RngIntElt
413{}
414   E := Curve(Q);
415   t := t_val(Q);
416   time sigma := formal_sigma_in_s2(E, prec);
418   print "s1 = ", s1;
419   time s2 := solve_for_s2(sigma, p);
420   print "s2 = ", s2;
421   time s2b := compute_s2(E, p, 40, s1);
422   print "s2b = ", s2b;
423   print "s2 - s2b = ", Parent(s2)!(s2 - s2b);
424   v := Evaluate(Evaluate(sigma, t),s2);
425   return v;
426end intrinsic;
427
428intrinsic height_of_point_modp(P::PtEll, p::RngIntElt) -> .
429{Returns p-adic height of point P.}
430   n := prep_multiple(P, p);
431   Q := n*P;
432   d := x_coord_d(Q);  // positive sqrt of denominator
433   HQ := extended_log(val_of_sigma(Q, p, 10) / d);
434   HP := HQ/n^2;
435   return GF(p)!(Integers()!HP);
436end intrinsic;
437
438intrinsic myExp(f::.) -> RngSerLaurElt
439{}
440   return &+[f^n/Factorial(n) : n in [0..AbsolutePrecision(f)]];
441end intrinsic;
442
443intrinsic myIsZero(a::.) -> .
444{}
445   return Valuation(a) ge Precision(Parent(a));
446end intrinsic;
447
448intrinsic myIntegral(f::RngSerLaurElt) -> RngSerLaurElt
449{
450 Integrate Laurent series, without stupid check on coefficient
451 of t^(-1) being 0, which doesn't work in MAGMA, since nonzero
452 for p-adic elements is messed up.
453}
454   require myIsZero(Coefficient(f,-1)) : "Coefficient of 1/t must be zero.";
455   pr := AbsolutePrecision(f);
456   t := Parent(f).1;
457   return &+[(Coefficient(f,n)/(n+1))*t^(n+1) : n in [0..pr-1] |
458                not myIsZero(Coefficient(f,n)) ] + O(t^pr);
459
460//   return -Coefficient(f,-2)*t^(-1) +
461//        &+[(Coefficient(f,n)/(n+1))*t^(n+1) : n in [0..pr-1]] + O(t^pr);
462end intrinsic;
463
464intrinsic sigma_alg2(E::CrvEll, p::RngIntElt, prec::RngIntElt :
465                    e2prec := 50) -> RngElt
466{
467   Compute and return the formal series for sigma.
468This doesn't work at all, because it requires integrating 1/t,
469which is impossible.
470}
471   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
473
474   // Compute w such that omega = w(t) dt
475   w := R!formal_omega(E, prec);
476   //t := Parent(w).1;
477
478   // Compute E2(E,omega)
479   e2 := E2(E, p, e2prec);
480
481   // Compute canonical X-function
482   x := R!formal_x(E, prec);
483   X := x + R!(a1^2 + 4*a2 + e2);
484
485   /* Solve the differential equation
486         -d/omega ( (d(sigma)/omega) * (1/sigma) ) = X
487         subject to omega = t + higher terms
488    */
489   intwX := myIntegral(w*X);
490print "w = ", w;
491print "intwX = ", intwX;
492print "w*intwX = ", w*intwX;
493   f := -myIntegral(w*intwX);
494print "f after int = ", f;
495   // But the constant term and linear term are not well-defined just from integrating,
496   // so we set them to 0, 1, respectively.
497   pr := AbsolutePrecision(f);
498   f := R!([0,1] cat [Coefficient(f,n) : n in [2..pr-1]]) + O(t^(pr-5));
499   print "f = ", f;
500/*   S<t> := PowerSeriesRing(BaseRing(R));
501   pr := AbsolutePrecision(f);
502Eltseq(f);
503print pr;
504   f := S!([0] cat [Coefficient(f,n) : n in [1..pr-1]]) + O(t^pr);
505*/
506   // Finally find sigma.
507print "f = ", f;
508   sigma := myExp(f);
509
510   return sigma;
511end intrinsic;
512
513intrinsic sigma_alg1(E::CrvEll, p::RngIntElt, prec::RngIntElt) -> .
514{}
515   time sigma_in_s2 := formal_sigma_in_s2(E, prec);
517   time s2 := solve_for_s2(sigma_in_s2, p);
519   pr := AbsolutePrecision(sigma_in_s2);
520   sigma := &+[Evaluate(Coefficient(sigma_in_s2,n),s2)*t^n : n in [1..pr-1]];
521   return sigma;
522end intrinsic;
523
524intrinsic sigma_using_e2(E::CrvEll, p::RngIntElt, prec::RngIntElt :
525                    e2prec := 50) -> .
526{
527   Compute and return the formal series for sigma.
528}
529   prec +:= 2;
530
531   // 1. Compute value of p-adic E_2, using Kedlaya's MW algorithm
532   e2 := E2(E,p,e2prec);
533
534   // 2. Define some notation.
536   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
537   S<t> := LaurentSeriesRing(K);
538
539   // 3. Compute formal Weierstrass p function in terms of t.
540   x  := S!formal_x(E, prec);
541   wp := x + (a1^2 + 4*a2)/12;
542
543   // 4. Compute wp(z), where z = FormalLog(t)
544   z_of_t := S!formal_log(E, prec);
545   t_of_z := Reverse(z_of_t);
546   wp_of_z := Evaluate(wp, t_of_z);
547   R<z> := S;
548
549   // 5. Compute
550   pr := AbsolutePrecision(wp_of_z);
551   // Let g be    1/z^2 - \wp + E2/12.  Notice that we compute this
552   // by just ignoring the coefficients of wp_of_z for n=-2.
553   g := e2/12 +  &+[z^n*Coefficient(-wp_of_z,n) : n in [-1..pr-1]] + O(z^pr);
554   /* It's *really* weird that it isn't -E2/12, since the canonical
555      eks function is \wp + E2/12.
556
557      There's a factor of -1 in my definition of E2 in the other file,
558      which I think Tate and I got out of Katz's paper.   Maybe that
559      sign is wrong?   There's some -1 normalization that is different
560      in two papers...
561   */
562
563   // 6. Integrate twice, exponentiate, etc.
564   sigma_of_z := z*myExp(myIntegral(myIntegral(g)));
565   sigma_of_t := Evaluate(sigma_of_z, z_of_t);
566   R<t> := PowerSeriesRing(K);
567
568   return R!sigma_of_t + O(t^(prec-2));
569end intrinsic;
570
571intrinsic height_function(E::CrvEll, p::RngIntElt, prec::RngIntElt) -> .
572{Returns p-adic height of point P.}
573   sig := sigma_using_e2(E, p, prec : e2prec := prec) ;
574   function h(P)
575      assert Curve(P) eq E;
576      n := prep_multiple(P, p);
577      Q := n*P;
578      d := x_coord_d(Q);  // positive sqrt of denominator
579      t := t_val(Q);
580      v := Evaluate(sig,t);
581      HQ := extended_log(v/d);
582      HP := HQ/n^2;
583      return HP/p;
584   end function;
585   return h;
586end intrinsic;
587
588intrinsic regulator(E::CrvEll, p::RngIntElt) -> .
589{}
590    return regulator(E,p,20);
591end intrinsic;
592
593intrinsic regulator(E::CrvEll, p::RngIntElt, prec::RngIntElt) -> .
594{
595Computes and returns the p-adic regulator of E.
596}
597   h := height_function(E, p, prec);
598   G, f := MordellWeilGroup(E);
599   B := [f(G.i) : i in [1..Ngens(G)] | Order(G.i) eq 0];
600   function pairing(P,Q)
601      return (h(P+Q)-h(P)-h(Q))/2;
602   end function;
603   K := Parent(h(B[1]));
604   M := MatrixAlgebra(K,#B)!0;
605   for i in [1..#B] do
606      for j in [i..#B] do
607         aij := pairing(B[i],B[j]);
608         M[i,j] := aij;
609         M[j,i] := aij;
610      end for;
611   end for;
612   return Determinant(M), M;
613end intrinsic;
614
615
616intrinsic val_of_sigma_alg2(Q::PtEll, p::RngIntElt,
617                       prec::RngIntElt) -> RngElt
618{
619  Return the value of the p-adic sigma function on Q, using the new algorithm.
620}
621   E := Curve(Q);
622   sigma := sigma_alg2(E, p, prec);
623   t := t_val(Q);
624   return Evaluate(sigma, t);
625end intrinsic;
626
627intrinsic height_of_point(P::PtEll, p::RngIntElt, prec::RngIntElt) -> .
628{Returns p-adic height of point P.}
629   n := prep_multiple(P, p);
630   Q := n*P;
631   d := x_coord_d(Q);  // positive sqrt of denominator
632   v := val_of_sigma(Q, p, prec);
633   HQ := extended_log(v/d);
634   HP := HQ/n^2;
635   return HP/p;
636end intrinsic;
637
638intrinsic experiment(p::RngIntElt, a_range::SeqEnum, prec::RngIntElt)
639{}
640   for a in a_range do
641      P, d, t := point(a);
642      if not t or d mod p eq 0 then
643         continue;
644      end if;
645      h := height_of_point(P, p, prec) ;
646      E := Curve(P);
647      printf "[%o, %o], ", a, h;
648   end for;
649end intrinsic;
650
651intrinsic ranks(p::RngIntElt, a_range::SeqEnum, prec::RngIntElt)
652{}
653   for a in a_range do
654      P, d, t := point(a);
655      if not t or d mod p eq 0 then
656         continue;
657      end if;
658      E := Curve(P);
659print aInvariants(E);
660      printf "a = %o\t\trank = %o\n", a, Rank(E);
661   end for;
662end intrinsic;
663
665{}
666   z := Integers()!x;
667   p := Prime(Parent(x));
668   i := 0;
669   s := "";
670   while z ne 0 and i lt Precision(Parent(x)) do
671      c := z mod p;
672      if c ne 0 then
673         if c ne 1 then
674            s *:= Sprintf("%o*", c);
675         end if;
676         s *:= Sprintf("%o^%o + ", p, i);
677      end if;
678      z := z div p;
679      i := i + 1;
680   end while;
681   s *:= Sprintf("O(%o^%o)", p, i);
682   return s;
683end intrinsic;
684
685
686