/******************************************************************* ******************************************************************* KEDLAYA'S ALGORITHM (in ODD CHARACTERISTIC) Mike Harrison. ******************************************************************* *******************************************************************/ LINEAR_LIMIT := 7; ALWAYS_CUBE := false; function GetRing(R1,Q,prec) L := LaurentSeriesRing(ChangePrecision(R1,prec)); P1 := PolynomialRing(L); return quo; end function; function myXGCD(p1,p2) // p1 and p2 are relatively prime integral polynomials over an // unramified pAdic field of bounded precision n. // It is assumed that the leading coeff of p1 is prime to p. // The function calculates and returns integral polynomials of // the same precision s.t. A*p1+B*p2 = 1 mod p^n. // Begins by lifting the Xgcd result over the Residue field and // iterates up mod p^i for 1 <= i <= n. R := Parent(p1); F,mp := ResidueClassField(IntegerRing(BaseRing(R))); S := PolynomialRing(F); p1r := S![mp(j) : j in Coefficients(p1)]; p2r := S![mp(j) : j in Coefficients(p2)]; _,Ar,Br := Xgcd(p1r,p2r); u := R!Ar; v := R!Br; A := u; B := v; delta := R!-1; p := Characteristic(F); for i in [1..BaseRing(R)`DefaultPrecision-1] do delta := (delta+u*p1+v*p2)/p; delr := S![mp(j) : j in Coefficients(delta)]; v1 := (-Br*delr) mod p1r; v := R!v1; u := R!((-(delr+v1*p2r)) div p1r); /* NB. the simpler u := R!((-Ar*delr) mod p2r); v := R!((-Br*delr) mod p1r); doesn't necessarily work if p2 has leading coeff div by p, when deg p2r < deg p2. In this case, u*p1+v*p2 != -delta mod p if deg p1r+deg p2r <= deg delr (< deg p1+deg p2) */ A +:= u*(p^i); B +:= v*(p^i); end for; return A,B; end function; function SigmaPowMat(M,m,n) // returns s^(n-1)(M)*s^(n-2)(M)*..*s(M)*M where M is // an m*m matrix over an unramified pAdic field and // s is the absolute Frobenius of that field. n >= 1. // Uses a basic left-right binary powering algorithm. if n eq 1 then return M; end if; //special case bits := Intseq(n,2); N := M; r := 1; for i := (#bits-1) to 1 by -1 do N := Matrix(m,[GaloisImage(x,r) : x in Eltseq(N)])*N; r *:= 2; if bits[i] eq 1 then N := Matrix(m,[GaloisImage(x,1) : x in Eltseq(N)])*M; r +:= 1; end if; end for; return N; end function; function FrobYInv(Q,p,N,x_pow,S,cube) // Computes p*(Frob y)^-1 (cube=false) or p*(Frob y)^(-3) (cube=true) // mod p^N. // This means calculating // (1+((FrobQ)(x^p)-(Q(x))^p)*X^p)^(-1/2){resp(-3/2)} // to the required precision in S (S defined in Kedlaya fn). // // Starts by computing terms in the binomial expansion of the above // then uses Newton iteration. The first part computes the most // appropriate number <= LINEAR_LIMIT of terms for the Newton phase. // In the Newton phase, the power series (in T) coefficients of powers // of x are truncated noting that result mod p^t contains no non-zero // coefficients of T beyond T^(p*(t-1)). // // Q(x) is the defining polynomial for the hyperelliptic curve (y^2=Q(x)) // and x_pow = x^(p-1) in S. R1 := BaseRing(BaseRing(S)); T := BaseRing(S).1; E := 1 + (T^p)*(Evaluate(Parent(Q1)![GaloisImage(j,1): j in Coefficients(Q1)], x_pow*S.1) - Evaluate(Q1,S.1)^p) where Q1 is PolynomialRing(R1)!Q; // now compute E^(-1/2) (E^(-3/2) if cube) mod p^N // ( this is weaker than mod T^(p(N-1)+1) ) // by "linear lift" followed by Newton iteration. //first backtrack to find starting point for Newton phase prec := N; adjs := [Booleans()|]; while prec gt LINEAR_LIMIT do Append (~adjs,IsOdd(prec)); if p eq 3 then prec := prec div 2; else prec := (prec+1) div 2; end if; end while; // now do the linear part Sc := GetRing(R1,Q,prec); Rc := BaseRing(BaseRing(Sc)); u := Sc!1; Epow := u; E1 := Sc![BaseRing(Sc)!a : a in Coefficients(E)]; half := cube select -(1/2) else (1/2); bincoeff := 1/1; for i in [1..prec] do bincoeff *:= (half-i)/i; Epow := (E1-1)*Epow; u := u + (Rc!bincoeff)*Epow; end for; delete Epow; // u is now the answer mod p^prec. Finish by Newton iteration // x -> (3*x - E*x^3)/2. if cube then E := E^3; end if; half := BaseRing(Parent(Q))!(1/2); if prec eq N then PolR := PolynomialRing(BaseRing(S)); end if; while prec lt N do tyme := Cputime(); prec *:= 2; if adjs[#adjs] then prec +:= ((p eq 3) select 1 else -1); end if; Prune(~adjs); if prec lt N then Sc := GetRing(R1,Q,prec); E1 := Sc![BaseRing(Sc)!a : a in Coefficients(E)]; else Sc := S; E1 := E; end if; T := BaseRing(Sc).1; PolR := PolynomialRing(BaseRing(Sc)); u := Sc![BaseRing(Sc)!a : a in Coefficients(u)]; v := Sc![j+O(T^(p*prec-(p-1))) : j in Coefficients(PolR!(u^2))]; u := (3*u - E1*(u*v))*(BaseRing(BaseRing(Sc))!(1/2)); // remove O(T^r) terms u := Sc![ elt where coeffs,v is Coefficients(j) : j in Coefficients(PolR!u)]; vprintf JacHypCnt, 3: "Reached precision %o time = %o\n",prec,Cputime(tyme); end while; // now "clean" the result mod T^(pN-(p-1)) u := S![ elt where coeffs,v is Coefficients(j+O(T^((p*(N-1))+1))) : j in Coefficients(PolR!(p*u))]; return (u * T^(((cube select 3 else 1)*p) div 2)); end function; function Convert(powTx,bdu,bdl,K) // Takes a differential powTx*(dx/y) where powTx is of the form // p0(T) + p1(T)*x +... pr(T)*x^r // (T := 1/y^2, pi(T) are finite-tailed Laurent series) // and returns [Ar,A(r+1),...],r where Ai = Coefficients(ai) // powTx = ar(x)*T^r + a(r+1)(x)*T^(r+1) + ... (ai in K[x]). // K is the unramified pAdic coefficient field. // bdu,bdl are upper and lower bounds for non-zero powers of // T in {p0,p1,...}. vec := [PowerSequence(K)|]; found_start := false; start_adj := 0; for i in [bdl..bdu] do v1 := [K!Coefficient(ser,i) : ser in Coefficients(powTx)]; if not found_start then if &and[RelativePrecision(k) eq 0 : k in v1] then start_adj +:= 1; continue; else found_start := true; end if; end if; Append(~vec,v1); end for; return vec,bdl+start_adj; end function; function PrecDiv(pol,d) // Used by ReduceB to avoid losing padic precision when // dividing polynomial pol by integer d (p may divide d). // If d isn't a padic unit, arbitrary "noise" is added to // restore full precision after the division. K := BaseRing(Parent(pol)); pold := pol/d; if Valuation(d,Prime(K)) ne 0 then pold := Parent(pol)! [(RelativePrecision(x) eq 0) select O(UniformizingElement(K)^K`DefaultPrecision) else ChangePrecision(x,Max(0,K`DefaultPrecision-Valuation(x))) : x in Coefficients(pold)]; end if; return pold; end function; function ReduceA(polys,precs,n0) // Performs the reduction of // {a_(n0-1)(x)*y^(2*(n0-1)) + .. + a1(x)*y^2 + a0(x)}*(dx/y) // to canonical return_val *(dx/y) form. PolR := Parent(precs[1]); if IsEmpty(polys) then return PolR!0; end if; d := Degree(precs[1]); pol := PolR!polys[1]; N := BaseRing(PolR)`DefaultPrecision; for k := (n0-1) to 0 by -1 do pol := ((k le (n0-#polys)) select PolR!0 else PolR!(polys[n0+1-k])) + pol*precs[1]; for j := (2*d-1) to d by -1 do lc := Coefficient(pol,j); if RelativePrecision(lc) ne 0 then pol := pol - ((PolR.1)^(j-d))* (ChangePrecision(lc,N)/((2*k-1)*d+2*(j+1))) * ((2*k+1)*precs[2]*PolR.1+2*(j+1-d)*precs[1]); end if; end for; pol := PolR![Coefficient(pol,i):i in [0..(d-1)]]; end for; lc := Coefficient(pol,d-1); if RelativePrecision(lc) ne 0 then pol := pol - (ChangePrecision(lc,N)/d)*precs[2]; end if; return pol; end function; function ReduceB(polys,precs,n0,cube) // Performs the reduction of // {a1(x)*(1/y^2) + .. a_n0(x)*(1/y^2n0)}*(dx/y^r) // to canonical return_val *(dx/y^r) form. // r = 1 if cube = false, else r = 3. PolR := Parent(precs[1]); if IsEmpty(polys) then return PolR!0; end if; divcase := (Valuation(LeadingCoefficient(precs[2])) gt 0); pol := PolR!polys[#polys]; for k := (n0-1+#polys) to (cube select 2 else 1) by -1 do p1 := (pol*precs[4]) mod precs[1]; if divcase then p2 := (pol-p1*precs[2]) div precs[1]; else p2 := (pol*precs[3]) mod precs[2]; end if; pol := ((k le n0) select PolR!0 else PolR!(polys[k-n0])) + p2 + PrecDiv(2*Derivative(p1),2*k-1); end for; return pol; end function; function UpperCoeffs(M,g,ppow,e_val) // Calcs the sequence of the upper coefficients (x^g,..x^(2g)) // of CharPoly(M) using the trace method. The magma intrinsic // could be used but this is slightly more efficient as only // Trace(M),Trace(M^2),..Trace(M^g) need be calculated rather // than Trace(M),..,Trace(M^(2g)). // If Nrows(M) = 2*g+1 then the extra eigenvalue of M, e_val, // is discarded. (e_val = q (1) if cube = false (true)). // The sequence [r(v)] is returned where, for a p-adic int v, // r(v) is the integer n s.t.|n|; Embed(BaseRing(Parent(poly)),ResidueClassField(R1)); S := GetRing(R1,poly,N); x := S.1; //R := LaurentSeriesRing(R1); //S1 := PolynomialRing(R); precs := [PolynomialRing(K)!poly]; //S := quo; Append(~precs,Derivative(precs[1])); A,B := myXGCD(precs[1],precs[2]); // A,B satisfy A*Q+B*Q'=1 where Q is the lift of poly to char 0. Append(~precs,A); Append(~precs,B); //Stage 1 - compute p*x^(p-1)*(y^Frob)^-1[3] vprintf JacHypCnt, 2: "Computing (y^Frobenius)^-%o ...\n",cube select 3 else 1; tyme :=Cputime(); x_pow := x^(p-1); difl := FrobYInv(PolynomialRing(R1)!poly,p,N,x_pow,S,cube)*x_pow; x_pow := x_pow*x; vprintf JacHypCnt, 2: "Expansion time: %o\n",Cputime(tyme); //Stage 2 - loop to calculate the rows of the Frobenius matrix vprint JacHypCnt, 2: "Reducing differentials ..."; R1 := quo; M := RMatrixSpace(R1,d,d)!0; i := 1; boundu := p*N+(p div 2)-1; S1 := PolynomialRing(BaseRing(S)); vtime JacHypCnt, 2: while true do boundl := (p div 2) - Floor((i*p-1)/(d+1)); polys,bot := Convert(S1!difl,boundu,boundl,K); diffr := ReduceA([polys[k] : k in [1..Min(1-bot,#polys)]],precs,-bot)+ ReduceB([polys[k] : k in [Max(2-bot,1)..#polys]],precs,Max(bot,1),cube); M[i] := RSpace(R1,d)![R1!Coefficient(diffr,k) : k in [0..(d-1)]]; if i eq d then break; end if; i +:= 1; difl := difl * x_pow; end while; print "M = ", M; print "CharacteristicPolynomial of M = ", CharacteristicPolynomial(M); vprint JacHypCnt, 2: "Computing Frobenius matrix by powering ..."; vtime JacHypCnt, 2: M := SigmaPowMat(M,d,n); // Now change the precision to N1+Val(p,g!). The Val(p.. is needed // to add random noise as the characteristic polynomial calculation // uses the Trace method which involves dividing by 1,2,..g for the // top terms (later terms aren't needed) with a corresponding loss // in precision. N2 := N1 + Valuation(Factorial(g),p); if N2 gt N then ChangePrecision(~K,N2); end if; M := Matrix(K,M); if N2 gt N1 then M := Matrix(K,d,[ChangePrecision(K!x,N2-Valuation(x)) : x in Eltseq(M)]); end if; tyme := Cputime(); q := #ResidueClassField(Integers(K)); UCoeffs := UpperCoeffs(M,g,p^N1,cube select 1 else q); CharP := PolynomialRing(IntegerRing())! ([UCoeffs[i]*q^(g+1-i): i in [1..g]] cat Reverse(UCoeffs)); vprintf JacHypCnt,3: "Characteristic polynomial time: %o\n",Cputime(tyme); return Parent(CharP)!Reverse(Coefficients(CharP)); end function; intrinsic kedlaya(poly::., p::., n::.) -> . {} return Kedlaya(poly, p, n); end intrinsic; function KedCharPolyOdd(C) Fld := BaseField(C); if Type(Fld) ne FldFin then error "The curve must be defined over a finite field!"; end if; p := Characteristic(Fld); if p eq 2 then error "Sorry! The kedlaya method can't currently handle char 2"; end if; n := Degree(Fld); C1 := SimplifiedModel(C); Q := HyperellipticPolynomials(C1); twist := false; lc := LeadingCoefficient(Q); if lc ne Fld!1 then if IsOdd(Degree(Q)) then Q := Evaluate(Q,lc*Parent(Q).1); Q := Q/(lc^(Degree(Q)+1)); else Q := Q/lc; if not IsSquare(lc) then twist := true; end if; end if; end if; vprint JacHypCnt, 1: "Applying Kedlaya's algorithm"; tyme := Cputime(); cpol := Kedlaya(Q,p,n); vprintf JacHypCnt, 1: "Total time: %o\n", Cputime(tyme); if twist then cpol := Evaluate(cpol,-(Parent(cpol).1)); end if; return cpol; end function; intrinsic kedlayacharpoly(C::.) -> . {} return KedCharPolyOdd(C); end intrinsic; /* Date: Thu, 8 Jul 2004 13:53:13 -0400 From: Nick Katz Subject: Re: convergence of the Eisenstein series of weight two To: mazur@math.harvard.edu, nmkatz@Math.Princeton.EDU Cc: tate@math.utexas.edu, was@math.harvard.edu It seems to me you want to use the interpretation of P as the "direction of the unit root subspace", that should make it fast to compute. Concretely, suppose we have a pair (E, \omega) over Z_p, and to fix ideas p is not 2 or 3. Then we write a Weierstrass eqn for E, y^2 = 4x^3 - g_2x - g_3, so that \omega is dx/y, and we denote by \eta the differential xdx/y. Then \omega and \eta form a Z_p basis of H^1_DR = H^1_cris, and the key step is to compute the matrix of absolute Frobenius (here Z_p linear, the advantage of working over Z_p: otherwise if over Witt vectors of an F_q, only \sigma-linear). [This calculation goes fast, because the matrix of Frobenius lives over the entire p-adic moduli space, and we are back in the glory days of Washnitzer- Monsky cohomology (of the open curve E - {origin}.] Okay, now suppose we have computed the matrix of Frob in the basis \omega, \eta. The unit root subspace is a direct factor, call it U, of the H^1, and we know that a complimentary direct factor is Fil^1 := the Z_p span of \omega. We also know that Frob(\omega) lies in pH^1, and this tells us that, mod p^n, U is the span of Frob^n(\eta). What this means concretely is that if we write, for each n, Frob^n(\eta) = a_n\omega + b_n\eta, then b_n is a unit (cong mod p to the n'th power of the hasse invariant) and that P is something like the ratio a_n/b_n (up to a sign and a factor 12 which i don't recall offhand but which is in my Antwerp appendix and also in my "p-adic interp. of real anal. Eis. series" paper. So in terms of speed of convergence, ONCE you have Frob, you have to iterate it n times to calculate P mod p^n. Best, Nick */ intrinsic padicprint(x::.) ->. {} z := Integers()!x; p := Prime(Parent(x)); i := 0; s := ""; while z ne 0 and i lt Precision(Parent(x)) do c := z mod p; if c ne 0 then if c ne 1 then s *:= Sprintf("%o*", c); end if; s *:= Sprintf("%o^%o + ", p, i); end if; z := z div p; i := i + 1; end while; s *:= Sprintf("O(%o^%o)", p, i); return s; end intrinsic; /* */ intrinsic TateCM_Examples() {} c4 := 2^9*11; c6 := 2^9*7*11^2; a4 := - c4/(2^4*3); a6 := - c6/(2^5*3^3); E := EllipticCurve([a4,a6]); print "sqrt(-11): mod5 ", (Integers()!E2(E,5,30)) mod 5; print "sqrt(-11): mod5^5 ", (Integers()!E2(E,5,30)) mod 5^5; c4 := 2^4*3*11; c6 := 2^6*3^3*11; a4 := - c4/(2^4*3); a6 := - c6/(2^5*3^3); E := EllipticCurve([a4,a6]); print "j=1: mod5 ", (Integers()!E2(E,5,30)) mod 5; print "j=1: mod5^5 ", (Integers()!E2(E,5,30)) mod 5^5; c4 := 2^5*11; c6 := -2^3*7*11^2; a4 := - c4/(2^4*3); a6 := - c6/(2^5*3^3); E := EllipticCurve([a4,a6]); print "j=2: mod5 ", (Integers()!E2(E,5,30)) mod 5; print "j=2: mod5^5 ", (Integers()!E2(E,5,30)) mod 5^5; c4 := -2^4*3; c6 := 0; a4 := - c4/(2^4*3); a6 := - c6/(2^5*3^3); E := EllipticCurve([a4,a6]); print "j=3: mod5 ", (Integers()!E2(E,5,30)) mod 5; print "j=3: mod5^5 ", (Integers()!E2(E,5,30)) mod 5^5; c4 := 2^5*3*19; c6 := -2^3*3^3*19^2; a4 := - c4/(2^4*3); a6 := - c6/(2^5*3^3); E := EllipticCurve([a4,a6]); print "j=4: mod5 ", (Integers()!E2(E,5,30)) mod 5; print "j=4: mod5^5 ", (Integers()!E2(E,5,30)) mod 5^5; end intrinsic; intrinsic Lagrange(points::SeqEnum) -> RngUPolElt { Input: points -- a sequence of length n>=1 of pairs (xi,yi=f(xi)), with the xi distinct Output: a polynomial f of degree <= n-1 that passes through those points. } n := #points; R := PolynomialRing(FieldOfFractions(Parent(points[1][1]))); function P(j) return points[j][2]* &*[(x-points[k][1])/(points[j][1]-points[k][1]) : k in [1..n] | k ne j]; end function; return &+[P(j) : j in [1..n]]; end intrinsic; intrinsic IsGoodOrdinary(E::CrvEll, p::RngIntElt) -> BoolElt {} return Conductor(E) mod p ne 0 and TraceOfFrobenius(ChangeRing(E,GF(p))) ne 0; end intrinsic; /* Family: y^2 = x^3 + x*t + t; over Q(t). */ intrinsic E2oft(t::RngIntElt, p::RngIntElt, prec::RngIntElt) -> FldPadElt {} E := EllipticCurve([t,t]); require IsGoodOrdinary(E,p) : "Fiber curve E_t must be good ordinary at argument 2."; return E2(E, p, prec); end intrinsic; intrinsic E2_def(p::RngIntElt, radius::RngIntElt, num_points::RngIntElt) -> . {} points := []; n := 1; while #points lt num_points do //t := p*Random(2,20) + p^radius*n; t := p^radius*n; n +:= Random(1,20); E := EllipticCurve([-1296+t,11664+t]); if IsGoodOrdinary(E,p) then e := E2(E,p,20); v := [t,e] cat [Integers()!e mod p^a : a in [1..radius+5]]; Append(~points, v); end if; end while; f := Lagrange(points); return f, points, [Valuation(c) : c in Eltseq(f)]; end intrinsic;