CoCalc Public Fileswww / talks / harvard-talk-2004-12-08 / e2heights / formal.mOpen with one click!
Author: William A. Stein
Compute Environment: Ubuntu 18.04 (Deprecated)
1
/*****************************************************************
2
3
4
5
6
****************************************************************/
7
8
Q := RationalField();
9
R<t> := LaurentSeriesRing(RationalField());
10
11
intrinsic 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;
24
end intrinsic;
25
26
intrinsic formal_x(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
27
{}
28
return t/formal_w(E,prec+5);
29
end intrinsic;
30
31
intrinsic formal_y(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
32
{}
33
return -1/formal_w(E,prec+6);
34
end intrinsic;
35
36
intrinsic formal_omega(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
37
{Compute the formal differential in terms of t, i.e., return the function f(t)
38
such that f(t)dt is the formal differerntial omega.
39
40
Use 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);
47
end intrinsic;
48
49
intrinsic 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
51
integral 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));
56
end intrinsic;
57
58
intrinsic 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);
64
end intrinsic;
65
66
intrinsic 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);
90
end intrinsic;
91
92
93
94
intrinsic 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
*/
143
end intrinsic;
144
145
intrinsic 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);
154
end intrinsic;
155
156
intrinsic 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));
169
end intrinsic;
170
171
intrinsic 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);
179
end 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
186
intrinsic 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)));
204
end intrinsic;
205
206
intrinsic 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;
221
end intrinsic;
222
223
224
225
/*******************************************************
226
227
Karl Rubin and I were wondering whether you could do a p-adic height
228
computation, the point of which is to figure out whether the p-adic
229
regulator of elliptic curves has a "tendency" to be not divisible by p. We
230
want to avoid a few degenerate cases, so the smallest p that is useful to
231
us is p=5. Is the strategy-- described below--, that systematically runs
232
through 5-adic height computations for points on quadratics twists of the
233
(non-3-Eisenstein) elliptic curve of conductor 37 reasonably easy to do?
234
235
236
%%%%%%%%%%%%%%%%%%%
237
Our untwisted equation is E: y^2 = 4x^3-4x+1. For any nonzero integer D we
238
consider the twisted elliptic curve E_D: D.y^2 = 4x^3-4x+1. Put F(x): =
239
4x^3-4x+1. Let x run through integers
240
241
a= 0, \pm 1, \pm 2, \pm3, ...
242
243
and for each of these a's,
244
245
1) Factor F(a) = d.b^2, where d is square-free;
246
247
2) throw out the ones where d is divisible by 5, or where the parity of the
248
rank of E_d is even;
249
250
3) for the rest, viewing P:=(a,b) as a rational point-- it will be of
251
infinite order--on the elliptic curve E_d, so one can compute the 5-adic
252
height of P on the elliptic curve E_d: this should be a 5-adic integer if
253
the height is normalized right, and there should be a good number of
254
instances where this height is a 5-adic unit; throw them out, and
255
256
4) print the rest.
257
%%%%%%%%%%%%%%%%
258
259
Note that we only want to know the 5-adic height modulo 5.
260
261
Then it would be our job to figure out whether strange things occur for
262
this printed list of instances: e.g., the point P might be divisible by 5,
263
or the Mordell-Weil rank might be > 2, or other stuff..., and we want to
264
know what other stuff there might be.
265
266
****************************************************/
267
268
intrinsic 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;
280
end intrinsic;
281
282
function prime_divisors(n)
283
return [F[1] : F in Factorization(n)];
284
end function;
285
286
intrinsic prep_multiple(P::PtEll, p::RngIntElt) -> PtEll, RngIntElt
287
{An integer n such that n*P is in the subgroup of
288
E that reduces to the identity mod p and lands in the connected
289
component 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;
297
end intrinsic;
298
299
intrinsic extended_log(x::FldPadElt) -> FldPadElt
300
{}
301
p := Prime(Parent(x));
302
u := x * p^(-Valuation(x));
303
return Log(u^(p-1))/(p-1);
304
end intrinsic;
305
306
intrinsic x_coord_d(Q::PtEll) -> RngIntElt
307
{}
308
t, d := IsSquare(Denominator(Q[1]));
309
assert t;
310
return d;
311
end intrinsic;
312
313
intrinsic t_val(Q::PtEll) -> FldRatElt
314
{Given input (x,y), return t = -x/y.}
315
return -Q[1]/Q[2];
316
end intrinsic;
317
318
intrinsic 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);
323
end intrinsic;
324
325
intrinsic 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;
332
end intrinsic;
333
334
intrinsic 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;
363
return pAdicField(p,n)!a;
364
end intrinsic;
365
366
intrinsic 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;
380
end intrinsic;
381
382
intrinsic 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);
388
s1 := pAdicField(p)!Coefficient(Coefficient(sigma,2),0);
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;
397
end intrinsic;
398
399
intrinsic 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);
407
end intrinsic;
408
409
intrinsic myExp(f::.) -> RngSerLaurElt
410
{}
411
return &+[f^n/Factorial(n) : n in [0..AbsolutePrecision(f)]];
412
end intrinsic;
413
414
intrinsic myIsZero(a::.) -> .
415
{}
416
return Valuation(a) ge Precision(Parent(a));
417
end intrinsic;
418
419
intrinsic 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);
433
end intrinsic;
434
435
intrinsic 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));
441
R<t> := LaurentSeriesRing(pAdicField(p));
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);
459
print "w = ", w;
460
print "intwX = ", intwX;
461
print "w*intwX = ", w*intwX;
462
f := -myIntegral(w*intwX);
463
print "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);
471
Eltseq(f);
472
print pr;
473
f := S!([0] cat [Coefficient(f,n) : n in [1..pr-1]]) + O(t^pr);
474
*/
475
// Finally find sigma.
476
print "f = ", f;
477
sigma := myExp(f);
478
479
return sigma;
480
end intrinsic;
481
482
intrinsic sigma_alg1(E::CrvEll, p::RngIntElt, prec::RngIntElt) -> .
483
{}
484
time sigma_in_s2 := formal_sigma_in_s2(E, prec);
485
s1 := pAdicField(p)!Coefficient(Coefficient(sigma_in_s2,2),0);
486
time s2 := solve_for_s2(sigma_in_s2, p);
487
R<t> := PowerSeriesRing(pAdicField(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;
491
end intrinsic;
492
493
intrinsic 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.
504
K := pAdicField(p,Precision(Parent(e2)));
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));
538
end intrinsic;
539
540
intrinsic 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;
555
end intrinsic;
556
557
intrinsic regulator(E::CrvEll, p::RngIntElt) -> .
558
{}
559
return regulator(E,p,20);
560
end intrinsic;
561
562
intrinsic regulator(E::CrvEll, p::RngIntElt, prec::RngIntElt) -> .
563
{
564
Computes 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;
583
end intrinsic;
584
585
intrinsic 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];
591
end intrinsic;
592
593
intrinsic 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);
602
end intrinsic;
603
604
intrinsic 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;
613
end intrinsic;
614
615
intrinsic 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;
626
end intrinsic;
627
628
intrinsic 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);
636
print aInvariants(E);
637
printf "a = %o\t\trank = %o\n", a, Rank(E);
638
end for;
639
end intrinsic;
640
641
intrinsic padicprint(x::.) -> .
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;
660
end intrinsic;
661
662
663