Sharedwww / talks / padic_heights / kedlaya.mOpen in CoCalc
Author: William A. Stein
1
/*******************************************************************
2
*******************************************************************
3
4
KEDLAYA'S ALGORITHM (in ODD CHARACTERISTIC)
5
6
Mike Harrison.
7
8
*******************************************************************
9
*******************************************************************/
10
11
LINEAR_LIMIT := 7;
12
ALWAYS_CUBE := false;
13
14
function GetRing(R1,Q,prec)
15
16
L<T> := LaurentSeriesRing(ChangePrecision(R1,prec));
17
P1 := PolynomialRing(L);
18
return quo<P1|P1!Q-T^-1>;
19
20
end function;
21
22
function myXGCD(p1,p2)
23
24
// p1 and p2 are relatively prime integral polynomials over an
25
// unramified pAdic field of bounded precision n.
26
// It is assumed that the leading coeff of p1 is prime to p.
27
// The function calculates and returns integral polynomials of
28
// the same precision s.t. A*p1+B*p2 = 1 mod p^n.
29
// Begins by lifting the Xgcd result over the Residue field and
30
// iterates up mod p^i for 1 <= i <= n.
31
R := Parent(p1);
32
F,mp := ResidueClassField(IntegerRing(BaseRing(R)));
33
S := PolynomialRing(F);
34
p1r := S![mp(j) : j in Coefficients(p1)];
35
p2r := S![mp(j) : j in Coefficients(p2)];
36
_,Ar,Br := Xgcd(p1r,p2r);
37
u := R!Ar; v := R!Br;
38
A := u; B := v;
39
delta := R!-1;
40
p := Characteristic(F);
41
for i in [1..BaseRing(R)`DefaultPrecision-1] do
42
delta := (delta+u*p1+v*p2)/p;
43
delr := S![mp(j) : j in Coefficients(delta)];
44
v1 := (-Br*delr) mod p1r;
45
v := R!v1;
46
u := R!((-(delr+v1*p2r)) div p1r);
47
/* NB. the simpler
48
u := R!((-Ar*delr) mod p2r);
49
v := R!((-Br*delr) mod p1r);
50
doesn't necessarily work if p2 has leading
51
coeff div by p, when deg p2r < deg p2.
52
In this case, u*p1+v*p2 != -delta mod p if
53
deg p1r+deg p2r <= deg delr (< deg p1+deg p2)
54
*/
55
A +:= u*(p^i);
56
B +:= v*(p^i);
57
end for;
58
return A,B;
59
60
end function;
61
62
function SigmaPowMat(M,m,n)
63
64
// returns s^(n-1)(M)*s^(n-2)(M)*..*s(M)*M where M is
65
// an m*m matrix over an unramified pAdic field and
66
// s is the absolute Frobenius of that field. n >= 1.
67
// Uses a basic left-right binary powering algorithm.
68
69
if n eq 1 then return M; end if; //special case
70
bits := Intseq(n,2);
71
N := M;
72
r := 1;
73
for i := (#bits-1) to 1 by -1 do
74
N := Matrix(m,[GaloisImage(x,r) : x in Eltseq(N)])*N;
75
r *:= 2;
76
if bits[i] eq 1 then
77
N := Matrix(m,[GaloisImage(x,1) : x in Eltseq(N)])*M;
78
r +:= 1;
79
end if;
80
end for;
81
return N;
82
83
end function;
84
85
function FrobYInv(Q,p,N,x_pow,S,cube)
86
87
// Computes p*(Frob y)^-1 (cube=false) or p*(Frob y)^(-3) (cube=true)
88
// mod p^N.
89
// This means calculating
90
// (1+((FrobQ)(x^p)-(Q(x))^p)*X^p)^(-1/2){resp(-3/2)}
91
// to the required precision in S (S defined in Kedlaya fn).
92
//
93
// Starts by computing terms in the binomial expansion of the above
94
// then uses Newton iteration. The first part computes the most
95
// appropriate number <= LINEAR_LIMIT of terms for the Newton phase.
96
// In the Newton phase, the power series (in T) coefficients of powers
97
// of x are truncated noting that result mod p^t contains no non-zero
98
// coefficients of T beyond T^(p*(t-1)).
99
//
100
// Q(x) is the defining polynomial for the hyperelliptic curve (y^2=Q(x))
101
// and x_pow = x^(p-1) in S.
102
103
R1 := BaseRing(BaseRing(S));
104
T := BaseRing(S).1;
105
E := 1 + (T^p)*(Evaluate(Parent(Q1)![GaloisImage(j,1):
106
j in Coefficients(Q1)], x_pow*S.1) - Evaluate(Q1,S.1)^p)
107
where Q1 is PolynomialRing(R1)!Q;
108
// now compute E^(-1/2) (E^(-3/2) if cube) mod p^N
109
// ( this is weaker than mod T^(p(N-1)+1) )
110
// by "linear lift" followed by Newton iteration.
111
112
//first backtrack to find starting point for Newton phase
113
prec := N;
114
adjs := [Booleans()|];
115
while prec gt LINEAR_LIMIT do
116
Append (~adjs,IsOdd(prec));
117
if p eq 3 then
118
prec := prec div 2;
119
else
120
prec := (prec+1) div 2;
121
end if;
122
end while;
123
// now do the linear part
124
Sc := GetRing(R1,Q,prec);
125
Rc := BaseRing(BaseRing(Sc));
126
u := Sc!1;
127
Epow := u;
128
E1 := Sc![BaseRing(Sc)!a : a in Coefficients(E)];
129
half := cube select -(1/2) else (1/2);
130
bincoeff := 1/1;
131
for i in [1..prec] do
132
bincoeff *:= (half-i)/i;
133
Epow := (E1-1)*Epow;
134
u := u + (Rc!bincoeff)*Epow;
135
end for;
136
delete Epow;
137
// u is now the answer mod p^prec. Finish by Newton iteration
138
// x -> (3*x - E*x^3)/2.
139
if cube then E := E^3; end if;
140
half := BaseRing(Parent(Q))!(1/2);
141
if prec eq N then PolR := PolynomialRing(BaseRing(S)); end if;
142
while prec lt N do
143
tyme := Cputime();
144
prec *:= 2;
145
if adjs[#adjs] then
146
prec +:= ((p eq 3) select 1 else -1);
147
end if;
148
Prune(~adjs);
149
if prec lt N then
150
Sc := GetRing(R1,Q,prec);
151
E1 := Sc![BaseRing(Sc)!a : a in Coefficients(E)];
152
else
153
Sc := S; E1 := E;
154
end if;
155
T := BaseRing(Sc).1;
156
PolR := PolynomialRing(BaseRing(Sc));
157
u := Sc![BaseRing(Sc)!a : a in Coefficients(u)];
158
v := Sc![j+O(T^(p*prec-(p-1))) : j in Coefficients(PolR!(u^2))];
159
u := (3*u - E1*(u*v))*(BaseRing(BaseRing(Sc))!(1/2));
160
// remove O(T^r) terms
161
u := Sc![ elt<Parent(j)|v,coeffs,-1> where coeffs,v is
162
Coefficients(j) : j in Coefficients(PolR!u)];
163
vprintf JacHypCnt, 3:
164
"Reached precision %o time = %o\n",prec,Cputime(tyme);
165
end while;
166
// now "clean" the result mod T^(pN-(p-1))
167
u := S![ elt<Parent(j)|v,coeffs,-1> where coeffs,v is
168
Coefficients(j+O(T^((p*(N-1))+1))) : j in Coefficients(PolR!(p*u))];
169
return (u * T^(((cube select 3 else 1)*p) div 2));
170
171
end function;
172
173
function Convert(powTx,bdu,bdl,K)
174
// Takes a differential powTx*(dx/y) where powTx is of the form
175
// p0(T) + p1(T)*x +... pr(T)*x^r
176
// (T := 1/y^2, pi(T) are finite-tailed Laurent series)
177
// and returns [Ar,A(r+1),...],r where Ai = Coefficients(ai)
178
// powTx = ar(x)*T^r + a(r+1)(x)*T^(r+1) + ... (ai in K[x]).
179
// K is the unramified pAdic coefficient field.
180
// bdu,bdl are upper and lower bounds for non-zero powers of
181
// T in {p0,p1,...}.
182
183
vec := [PowerSequence(K)|];
184
found_start := false;
185
start_adj := 0;
186
for i in [bdl..bdu] do
187
v1 := [K!Coefficient(ser,i) : ser in Coefficients(powTx)];
188
if not found_start then
189
if &and[RelativePrecision(k) eq 0 : k in v1] then
190
start_adj +:= 1;
191
continue;
192
else
193
found_start := true;
194
end if;
195
end if;
196
Append(~vec,v1);
197
end for;
198
return vec,bdl+start_adj;
199
200
end function;
201
202
function PrecDiv(pol,d)
203
204
// Used by ReduceB to avoid losing padic precision when
205
// dividing polynomial pol by integer d (p may divide d).
206
// If d isn't a padic unit, arbitrary "noise" is added to
207
// restore full precision after the division.
208
209
K := BaseRing(Parent(pol));
210
pold := pol/d;
211
if Valuation(d,Prime(K)) ne 0 then
212
pold := Parent(pol)!
213
[(RelativePrecision(x) eq 0) select
214
O(UniformizingElement(K)^K`DefaultPrecision) else
215
ChangePrecision(x,Max(0,K`DefaultPrecision-Valuation(x))) :
216
x in Coefficients(pold)];
217
end if;
218
return pold;
219
220
end function;
221
222
function ReduceA(polys,precs,n0)
223
224
// Performs the reduction of
225
// {a_(n0-1)(x)*y^(2*(n0-1)) + .. + a1(x)*y^2 + a0(x)}*(dx/y)
226
// to canonical return_val *(dx/y) form.
227
228
PolR := Parent(precs[1]);
229
if IsEmpty(polys) then
230
return PolR!0;
231
end if;
232
d := Degree(precs[1]);
233
pol := PolR!polys[1];
234
N := BaseRing(PolR)`DefaultPrecision;
235
for k := (n0-1) to 0 by -1 do
236
pol := ((k le (n0-#polys)) select PolR!0 else PolR!(polys[n0+1-k])) +
237
pol*precs[1];
238
for j := (2*d-1) to d by -1 do
239
lc := Coefficient(pol,j);
240
if RelativePrecision(lc) ne 0 then
241
pol := pol - ((PolR.1)^(j-d))*
242
(ChangePrecision(lc,N)/((2*k-1)*d+2*(j+1))) *
243
((2*k+1)*precs[2]*PolR.1+2*(j+1-d)*precs[1]);
244
end if;
245
end for;
246
pol := PolR![Coefficient(pol,i):i in [0..(d-1)]];
247
end for;
248
lc := Coefficient(pol,d-1);
249
if RelativePrecision(lc) ne 0 then
250
pol := pol - (ChangePrecision(lc,N)/d)*precs[2];
251
end if;
252
return pol;
253
254
end function;
255
256
function ReduceB(polys,precs,n0,cube)
257
258
// Performs the reduction of
259
// {a1(x)*(1/y^2) + .. a_n0(x)*(1/y^2n0)}*(dx/y^r)
260
// to canonical return_val *(dx/y^r) form.
261
// r = 1 if cube = false, else r = 3.
262
263
PolR := Parent(precs[1]);
264
if IsEmpty(polys) then
265
return PolR!0;
266
end if;
267
divcase := (Valuation(LeadingCoefficient(precs[2])) gt 0);
268
pol := PolR!polys[#polys];
269
for k := (n0-1+#polys) to (cube select 2 else 1) by -1 do
270
p1 := (pol*precs[4]) mod precs[1];
271
if divcase then
272
p2 := (pol-p1*precs[2]) div precs[1];
273
else
274
p2 := (pol*precs[3]) mod precs[2];
275
end if;
276
pol := ((k le n0) select PolR!0 else PolR!(polys[k-n0])) +
277
p2 + PrecDiv(2*Derivative(p1),2*k-1);
278
end for;
279
return pol;
280
281
end function;
282
283
function UpperCoeffs(M,g,ppow,e_val)
284
285
// Calcs the sequence of the upper coefficients (x^g,..x^(2g))
286
// of CharPoly(M) using the trace method. The magma intrinsic
287
// could be used but this is slightly more efficient as only
288
// Trace(M),Trace(M^2),..Trace(M^g) need be calculated rather
289
// than Trace(M),..,Trace(M^(2g)).
290
// If Nrows(M) = 2*g+1 then the extra eigenvalue of M, e_val,
291
// is discarded. (e_val = q (1) if cube = false (true)).
292
// The sequence [r(v)] is returned where, for a p-adic int v,
293
// r(v) is the integer n s.t.|n|<ppow/2 and n=v mod ppow.
294
295
pAd := pAdicField(Parent(M[1,1]));
296
N := M;
297
UCs := [pAd!1]; // coeffs (highest to lowest) of CharPoly(M)
298
TrPows := [pAd|]; // [Trace(M),Trace(M^2),...]
299
for i in [1..g] do
300
Append(~TrPows,Eltseq(Trace(N))[1]);
301
Append(~UCs, (- &+[TrPows[j]*UCs[i+1-j] : j in [1..i]])/i);
302
if i lt g then N := N*M; end if;
303
end for;
304
if Nrows(M) ne 2*g then // original Q(x) of even degree
305
for i in [1..g] do
306
UCs[i+1] := UCs[i+1]+e_val*UCs[i];
307
end for;
308
end if;
309
return [((2*uc gt ppow) select uc-ppow else uc) where uc is
310
(IntegerRing()!x) mod ppow : x in UCs];
311
312
end function;
313
314
/*
315
if p eq 0 or n eq 0 then
316
R := Parent(poly);
317
k := BaseRing(R);
318
p := Characteristic(k);
319
n := Degree(k);
320
end if;
321
*/
322
323
function Kedlaya(poly, p, n)
324
// Main function for the (odd char) Kedlaya algorithm.
325
// Computes the numerator of the zeta function for
326
// y^2 = poly (defined over over GF(p^n)).
327
// The cube boolean variable indicates which differential
328
// basis is used for cohomological calculations -
329
// (dx/y), x(dx/y), x^2(dx/y), ... if cube = false
330
// (dx/y^3), x(dx/y^3), ... if cube = true
331
// The 1st is the natural basis, but leads to a non-integral
332
// transformation matrix if p is small compared to the genus.
333
// The strategy is to use the first if this is not the case
334
// unless the ALWAYS_CUBE value is true
335
336
d := Degree(poly)-1;
337
g := d div 2;
338
cube := true;
339
if not ALWAYS_CUBE then
340
if (IsEven(d) and p ge d) or // deg=2*g+1
341
(IsOdd(d) and p gt g) then // deg=2*g+2
342
cube := false;
343
end if;
344
end if;
345
346
N1 := Ceiling((n*g/2)+Log(p,2*Binomial(2*g,g)));
347
N := N1 + Floor(Log(p,2*N1))+1;
348
K := UnramifiedExtension(pAdicField(p,N),n);
349
R1 := quo<Integers(K)|UniformizingElement(K)^N>;
350
Embed(BaseRing(Parent(poly)),ResidueClassField(R1));
351
S := GetRing(R1,poly,N);
352
x := S.1;
353
//R<T> := LaurentSeriesRing(R1);
354
//S1<x> := PolynomialRing(R);
355
precs := [PolynomialRing(K)!poly];
356
//S<x> := quo<S1|S1!poly-T^-1>;
357
Append(~precs,Derivative(precs[1]));
358
A,B := myXGCD(precs[1],precs[2]);
359
// A,B satisfy A*Q+B*Q'=1 where Q is the lift of poly to char 0.
360
Append(~precs,A);
361
Append(~precs,B);
362
363
//Stage 1 - compute p*x^(p-1)*(y^Frob)^-1[3]
364
vprintf JacHypCnt, 2:
365
"Computing (y^Frobenius)^-%o ...\n",cube select 3 else 1;
366
tyme :=Cputime();
367
x_pow := x^(p-1);
368
difl := FrobYInv(PolynomialRing(R1)!poly,p,N,x_pow,S,cube)*x_pow;
369
x_pow := x_pow*x;
370
vprintf JacHypCnt, 2: "Expansion time: %o\n",Cputime(tyme);
371
372
//Stage 2 - loop to calculate the rows of the Frobenius matrix
373
vprint JacHypCnt, 2: "Reducing differentials ...";
374
R1 := quo<Integers(K)|UniformizingElement(K)^N1>;
375
M := RMatrixSpace(R1,d,d)!0;
376
i := 1;
377
boundu := p*N+(p div 2)-1;
378
S1 := PolynomialRing(BaseRing(S));
379
vtime JacHypCnt, 2:
380
while true do
381
boundl := (p div 2) - Floor((i*p-1)/(d+1));
382
polys,bot := Convert(S1!difl,boundu,boundl,K);
383
diffr := ReduceA([polys[k] : k in [1..Min(1-bot,#polys)]],precs,-bot)+
384
ReduceB([polys[k] : k in [Max(2-bot,1)..#polys]],precs,Max(bot,1),cube);
385
M[i] := RSpace(R1,d)![R1!Coefficient(diffr,k) : k in [0..(d-1)]];
386
if i eq d then break; end if;
387
i +:= 1;
388
difl := difl * x_pow;
389
end while;
390
print "M = ", M;
391
print "CharacteristicPolynomial of M = ", CharacteristicPolynomial(M);
392
393
vprint JacHypCnt, 2: "Computing Frobenius matrix by powering ...";
394
vtime JacHypCnt, 2:
395
M := SigmaPowMat(M,d,n);
396
// Now change the precision to N1+Val(p,g!). The Val(p.. is needed
397
// to add random noise as the characteristic polynomial calculation
398
// uses the Trace method which involves dividing by 1,2,..g for the
399
// top terms (later terms aren't needed) with a corresponding loss
400
// in precision.
401
N2 := N1 + Valuation(Factorial(g),p);
402
if N2 gt N then ChangePrecision(~K,N2); end if;
403
M := Matrix(K,M);
404
if N2 gt N1 then
405
M := Matrix(K,d,[ChangePrecision(K!x,N2-Valuation(x)) : x in Eltseq(M)]);
406
end if;
407
tyme := Cputime();
408
q := #ResidueClassField(Integers(K));
409
UCoeffs := UpperCoeffs(M,g,p^N1,cube select 1 else q);
410
CharP := PolynomialRing(IntegerRing())!
411
([UCoeffs[i]*q^(g+1-i): i in [1..g]] cat Reverse(UCoeffs));
412
vprintf JacHypCnt,3: "Characteristic polynomial time: %o\n",Cputime(tyme);
413
return Parent(CharP)!Reverse(Coefficients(CharP));
414
415
end function;
416
417
intrinsic kedlaya(poly::., p::., n::.) -> .
418
{}
419
return Kedlaya(poly, p, n);
420
end intrinsic;
421
422
423
function KedCharPolyOdd(C)
424
425
Fld := BaseField(C);
426
if Type(Fld) ne FldFin then
427
error "The curve must be defined over a finite field!";
428
end if;
429
p := Characteristic(Fld);
430
if p eq 2 then
431
error "Sorry! The kedlaya method can't currently handle char 2";
432
end if;
433
n := Degree(Fld);
434
C1 := SimplifiedModel(C);
435
Q := HyperellipticPolynomials(C1);
436
twist := false;
437
lc := LeadingCoefficient(Q);
438
if lc ne Fld!1 then
439
if IsOdd(Degree(Q)) then
440
Q := Evaluate(Q,lc*Parent(Q).1);
441
Q := Q/(lc^(Degree(Q)+1));
442
else
443
Q := Q/lc;
444
if not IsSquare(lc) then twist := true; end if;
445
end if;
446
end if;
447
vprint JacHypCnt, 1: "Applying Kedlaya's algorithm";
448
tyme := Cputime();
449
cpol := Kedlaya(Q,p,n);
450
vprintf JacHypCnt, 1: "Total time: %o\n", Cputime(tyme);
451
if twist then
452
cpol := Evaluate(cpol,-(Parent(cpol).1));
453
end if;
454
return cpol;
455
456
end function;
457
458
intrinsic kedlayacharpoly(C::.) -> .
459
{}
460
return KedCharPolyOdd(C);
461
end intrinsic;
462
463
464
/*
465
Date: Thu, 8 Jul 2004 13:53:13 -0400
466
From: Nick Katz <nmk@Math.Princeton.EDU>
467
Subject: Re: convergence of the Eisenstein series of weight two
468
To: mazur@math.harvard.edu, nmkatz@Math.Princeton.EDU
469
Cc: tate@math.utexas.edu, was@math.harvard.edu
470
471
It seems to me you want to use the interpretation of P as the
472
"direction of the unit root subspace", that should make it fast to
473
compute. Concretely, suppose we have a pair (E, \omega) over Z_p, and
474
to fix ideas p is not 2 or 3. Then we write a Weierstrass eqn for E,
475
y^2 = 4x^3 - g_2x - g_3, so that \omega is dx/y, and we denote by \eta
476
the differential xdx/y. Then \omega and \eta form a Z_p basis of
477
H^1_DR = H^1_cris, and the key step is to compute the matrix of
478
absolute Frobenius (here Z_p linear, the advantage of working over
479
Z_p: otherwise if over Witt vectors of an F_q, only \sigma-linear).
480
[This calculation goes fast, because the matrix of Frobenius lives
481
over the entire p-adic moduli space, and we are back in the glory days
482
of Washnitzer- Monsky cohomology (of the open curve E - {origin}.]
483
484
Okay, now suppose we have computed the matrix of Frob in the
485
basis \omega, \eta. The unit root subspace is a direct factor, call it
486
U, of the H^1, and we know that a complimentary direct factor is Fil^1
487
:= the Z_p span of \omega. We also know that Frob(\omega) lies in
488
pH^1, and this tells us that, mod p^n, U is the span of
489
Frob^n(\eta). What this means concretely is that if we write, for each
490
n,
491
492
Frob^n(\eta) = a_n\omega + b_n\eta,
493
494
then b_n is a unit (cong mod p to the n'th power of the hasse
495
invariant) and that P is something like the ratio a_n/b_n (up to a
496
sign and a factor 12 which i don't recall offhand but which is in my
497
Antwerp appendix and also in my "p-adic interp. of real
498
anal. Eis. series" paper.
499
500
So in terms of speed of convergence, ONCE you have Frob, you
501
have to iterate it n times to calculate P mod p^n. Best, Nick
502
*/
503
504
intrinsic padicprint(x::.) ->.
505
{}
506
z := Integers()!x;
507
p := Prime(Parent(x));
508
i := 0;
509
s := "";
510
while z ne 0 and i lt Precision(Parent(x)) do
511
c := z mod p;
512
if c ne 0 then
513
if c ne 1 then
514
s *:= Sprintf("%o*", c);
515
end if;
516
s *:= Sprintf("%o^%o + ", p, i);
517
end if;
518
z := z div p;
519
i := i + 1;
520
end while;
521
s *:= Sprintf("O(%o^%o)", p, i);
522
return s;
523
end intrinsic;
524
525
526
527
528
/*
529
*/
530
531
intrinsic TateCM_Examples()
532
{}
533
c4 := 2^9*11;
534
c6 := 2^9*7*11^2;
535
a4 := - c4/(2^4*3);
536
a6 := - c6/(2^5*3^3);
537
E := EllipticCurve([a4,a6]);
538
print "sqrt(-11): mod5 ", (Integers()!E2(E,5,30)) mod 5;
539
print "sqrt(-11): mod5^5 ", (Integers()!E2(E,5,30)) mod 5^5;
540
541
c4 := 2^4*3*11;
542
c6 := 2^6*3^3*11;
543
a4 := - c4/(2^4*3);
544
a6 := - c6/(2^5*3^3);
545
E := EllipticCurve([a4,a6]);
546
print "j=1: mod5 ", (Integers()!E2(E,5,30)) mod 5;
547
print "j=1: mod5^5 ", (Integers()!E2(E,5,30)) mod 5^5;
548
549
c4 := 2^5*11;
550
c6 := -2^3*7*11^2;
551
a4 := - c4/(2^4*3);
552
a6 := - c6/(2^5*3^3);
553
E := EllipticCurve([a4,a6]);
554
print "j=2: mod5 ", (Integers()!E2(E,5,30)) mod 5;
555
print "j=2: mod5^5 ", (Integers()!E2(E,5,30)) mod 5^5;
556
557
c4 := -2^4*3;
558
c6 := 0;
559
a4 := - c4/(2^4*3);
560
a6 := - c6/(2^5*3^3);
561
E := EllipticCurve([a4,a6]);
562
print "j=3: mod5 ", (Integers()!E2(E,5,30)) mod 5;
563
print "j=3: mod5^5 ", (Integers()!E2(E,5,30)) mod 5^5;
564
565
c4 := 2^5*3*19;
566
c6 := -2^3*3^3*19^2;
567
a4 := - c4/(2^4*3);
568
a6 := - c6/(2^5*3^3);
569
E := EllipticCurve([a4,a6]);
570
print "j=4: mod5 ", (Integers()!E2(E,5,30)) mod 5;
571
print "j=4: mod5^5 ", (Integers()!E2(E,5,30)) mod 5^5;
572
573
end intrinsic;
574
575
intrinsic Lagrange(points::SeqEnum) -> RngUPolElt
576
{
577
Input:
578
points -- a sequence of length n>=1 of pairs (xi,yi=f(xi)), with the xi distinct
579
580
Output:
581
a polynomial f of degree <= n-1 that passes through those points.
582
}
583
n := #points;
584
R<x> := PolynomialRing(FieldOfFractions(Parent(points[1][1])));
585
function P(j)
586
return points[j][2]* &*[(x-points[k][1])/(points[j][1]-points[k][1]) : k in [1..n] | k ne j];
587
end function;
588
return &+[P(j) : j in [1..n]];
589
590
end intrinsic;
591
592
intrinsic IsGoodOrdinary(E::CrvEll, p::RngIntElt) -> BoolElt
593
{}
594
return Conductor(E) mod p ne 0 and
595
TraceOfFrobenius(ChangeRing(E,GF(p))) ne 0;
596
end intrinsic;
597
598
/*
599
Family:
600
601
y^2 = x^3 + x*t + t;
602
603
over Q(t).
604
*/
605
606
intrinsic E2oft(t::RngIntElt, p::RngIntElt, prec::RngIntElt) -> FldPadElt
607
{}
608
E := EllipticCurve([t,t]);
609
require IsGoodOrdinary(E,p) : "Fiber curve E_t must be good ordinary at argument 2.";
610
return E2(E, p, prec);
611
end intrinsic;
612
613
intrinsic E2_def(p::RngIntElt, radius::RngIntElt, num_points::RngIntElt) -> .
614
{}
615
points := [];
616
n := 1;
617
while #points lt num_points do
618
//t := p*Random(2,20) + p^radius*n;
619
t := p^radius*n;
620
n +:= Random(1,20);
621
E := EllipticCurve([-1296+t,11664+t]);
622
if IsGoodOrdinary(E,p) then
623
e := E2(E,p,20);
624
v := [t,e] cat [Integers()!e mod p^a : a in [1..radius+5]];
625
Append(~points, v);
626
end if;
627
end while;
628
f := Lagrange(points);
629
return f, points, [Valuation(c) : c in Eltseq(f)];
630
end intrinsic;
631