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