Author: William A. Stein
1/*******************************************************************
2 *******************************************************************
3
4      KEDLAYA'S ALGORITHM (in ODD CHARACTERISTIC)
5
6              Mike Harrison.
7
8*******************************************************************
9*******************************************************************/
10
11LINEAR_LIMIT := 7;
12ALWAYS_CUBE := false;
13
14function 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
20end function;
21
22function 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
60end function;
61
62function 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
83end function;
84
85function 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;
115    while prec gt LINEAR_LIMIT do
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;
146            prec +:= ((p eq 3) select 1 else -1);
147        end if;
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
171end function;
172
173function 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;
186    for i in [bdl..bdu] do
187        v1 := [K!Coefficient(ser,i) : ser in Coefficients(powTx)];
189	    if &and[RelativePrecision(k) eq 0 : k in v1] then
191		continue;
192	    else
193	        found_start := true;
194	    end if;
195	end if;
196	Append(~vec,v1);
197    end for;
199
200end function;
201
202function 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)^KDefaultPrecision) else
215           ChangePrecision(x,Max(0,KDefaultPrecision-Valuation(x))) :
216            x in Coefficients(pold)];
217    end if;
218    return pold;
219
220end function;
221
222function 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);
229    if IsEmpty(polys) then
230        return PolR!0;
231    end if;
232    d := Degree(precs);
233    pol := PolR!polys;
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;
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*PolR.1+2*(j+1-d)*precs);
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;
251    end if;
252    return pol;
253
254end function;
255
256function 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);
264    if IsEmpty(polys) then
265        return PolR!0;
266    end if;
267    divcase := (Valuation(LeadingCoefficient(precs)) 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) mod precs;
271        if divcase then
272           p2 := (pol-p1*precs) div precs;
273        else
274           p2 := (pol*precs) mod precs;
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
281end function;
282
283function 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
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)));
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
312end 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
323function 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;
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));
358    A,B := myXGCD(precs,precs);
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
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
415end function;
416
417intrinsic kedlaya(poly::., p::., n::.) -> .
418{}
419   return Kedlaya(poly, p, n);
420end intrinsic;
421
422
423function 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;
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
456end function;
457
458intrinsic kedlayacharpoly(C::.) -> .
459{}
460  return  KedCharPolyOdd(C);
461end intrinsic;
462
463
464/*
465Date: Thu, 8 Jul 2004 13:53:13 -0400
466From: Nick Katz <nmk@Math.Princeton.EDU>
467Subject: Re: convergence of the Eisenstein series of weight two
468To: mazur@math.harvard.edu, nmkatz@Math.Princeton.EDU
469Cc: tate@math.utexas.edu, was@math.harvard.edu
470
471It 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
473compute. Concretely, suppose we have a pair (E, \omega) over Z_p, and
474to fix ideas p is not 2 or 3.  Then we write a Weierstrass eqn for E,
475y^2 = 4x^3 - g_2x - g_3, so that \omega is dx/y, and we denote by \eta
476the differential xdx/y. Then \omega and \eta form a Z_p basis of
477H^1_DR = H^1_cris, and the key step is to compute the matrix of
478absolute Frobenius (here Z_p linear, the advantage of working over
479Z_p: otherwise if over Witt vectors of an F_q, only \sigma-linear).
480[This calculation goes fast, because the matrix of Frobenius lives
481over the entire p-adic moduli space, and we are back in the glory days
482of Washnitzer- Monsky cohomology (of the open curve E - {origin}.]
483
484        Okay, now suppose we have computed the matrix of Frob in the
485basis \omega, \eta. The unit root subspace is a direct factor, call it
486U, 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
488pH^1, and this tells us that, mod p^n, U is the span of
489Frob^n(\eta). What this means concretely is that if we write, for each
490n,
491
492   Frob^n(\eta) = a_n\omega + b_n\eta,
493
494then b_n is a unit (cong mod p to the n'th power of the hasse
495invariant) and that P is something like the ratio a_n/b_n (up to a
496sign and a factor 12 which i don't recall offhand but which is in my
497Antwerp appendix and also in my "p-adic interp. of real
498anal. Eis. series" paper.
499
500        So in terms of speed of convergence, ONCE you have Frob, you
501have to iterate it n times to calculate P mod p^n. Best, Nick
502*/
503
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;
523end intrinsic;
524
525
526
527
528/*
529*/
530
531intrinsic 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
573end intrinsic;
574
575intrinsic Lagrange(points::SeqEnum) -> RngUPolElt
576{
577Input:
578   points -- a sequence of length n>=1 of pairs (xi,yi=f(xi)), with the xi distinct
579
580Output:
581   a polynomial f of degree <= n-1 that passes through those points.
582}
583   n := #points;
584   R<x> := PolynomialRing(FieldOfFractions(Parent(points)));
585   function P(j)
586      return points[j]* &*[(x-points[k])/(points[j]-points[k]) : k in [1..n] | k ne j];
587   end function;
588   return &+[P(j) : j in [1..n]];
589
590end intrinsic;
591
592intrinsic IsGoodOrdinary(E::CrvEll, p::RngIntElt) -> BoolElt
593{}
594   return Conductor(E) mod p ne 0 and
595           TraceOfFrobenius(ChangeRing(E,GF(p))) ne 0;
596end intrinsic;
597
598/*
599Family:
600
601 y^2 = x^3 + x*t + t;
602
603over Q(t).
604*/
605
606intrinsic 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);
611end intrinsic;
612
613intrinsic 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;