\\ (public version of hell.gp) \\======================================================== \\ \\ p-adic height calculations \\ \\ This pari-gp script computes the canonical p-adic height \\ for an elliptic curve with good and ordinary reduction at p \\ \\ Author : chris wuthrich -> c.wuthrich@dpmms.cam.ac.uk \\ \\ Version : 25/5/2004 \\ \\ \\ usage : load the functions with \r hell.gp \\ \\ (!) this will kill some variables t w z a (that I need to have in a certain order) \\ \\ type ?function for more information \\ \\ hellinit(el,bigO) \\ computes the necessary information with precision bigO, \\ for having two significant digits for the heights one need about bigO = p ! \\ but even bigO=5 will give the first digit \\ \\ hellinitp(el,p,bigO) \\ computes the specific information at p \\ \\ hellatp(he,p) \\ as hellinitp, but starting from the output he from hellinit \\ \\ hellheight(hep,P) \\ computes the p-adic height of P. hep is the output of either hellinitp or hellatp \\ \\ hellbil(hep,P,Q) \\ computes the p-adic height pairing of P and Q (not optimized yet) \\ \\ hellreg(hep, ba) \\ computes the p-adic regulator with respect to the basis ba. \\ \\ helleuler(el, ba, p, flag=1) \\ prints information on the euler-characteristic of e for the basis ba \\ \\ example : \\ \r hell.gp \\ he=hellinit([0,0,1,-1,0],30); \\ hep=hellatp(he,5); \\ P=[0,0]; \\ hellheight(hep,P) \\ -> 3*5 + 3*5^3 + O(5^4) \\ hep=hellatp(he,13); \\ hellheight(hep,P) \\ -> 6*13^2 + 9*13^3 + O(13^4) \\ \\ e=ellinit([0,0,1,-1,0]); \\ helleuler(e,[[0,0]],13) \\ -> \\ curve of conductor 37 \\ reduction at p good \\ global torsion 13^0 \\ tamagawa 13^0 \\ regulator 13^2 \\ \\ Euler char/Sha 13^1 \\========================================================= \\ some variables need to appear in the right order if(hellisinitalized == 1, , print("First initalisation of hell.gp. I will kill the variables t,w,a and z"); kill(t); kill(w); kill(a); kill(z); w;t;a;); hellisinitalized=1; \\========================================================= \\ basic definitions for working with elliptic curves hellee(x)=if(length(x)==2,denominator(x[2])/denominator(x[1]),0) hellaa(x)=numerator(x[1]) hellbb(x)=numerator(x[2]) helltt(x)= if(x[2]==0,1,-x[1]/x[2]) \\ ellnp is the order of the group of non-singular points on the \\ reduction of el mod p hellnp(el,p)= { local(elll); elll=ellchangecurve(el,ellglobalred(el)[2]); return(if(el.disc % p,p+1-ellap(elll,p), p-ellap(elll,p)))} hellisgoodordinary(el,p)=if(el.disc % p && (p+1-hellnp(el,p)) % p,1,0); hellisanomalous(el,p)=(Mod(e.disc, p) != 0) && (Mod(hellnp(el,p), p)==0); hellbadprimes(el)=if(factor(el.disc)[1,1]==-1,vecextract(factor(el.disc)[,1]~,"^1"),factor(el.disc)[,1]~); \\gives back 1 if P reduces to a smooth point mod p hellissmooth(el,P,p) = { local(xp,yp,ep); if (el.disc %p, return(1), \\good reduction if(Mod(hellee(P),p)==0, return(1), \\zero section is always smooth ep=Mod(hellee(P),p); yp=Mod(hellbb(P),p)*ep^(-3); xp=Mod(hellaa(P),p)*ep^(-2); if((el.a1*yp == 3*xp^2+2*el.a2*xp+el.a4)&&(2*yp + el.a1*xp+el.a3 ==0), return(0), return(1)))) } \\ gives back 1 if P reduces to a smooth point mod p for all primes p hellisallsmooth(el,P)= { local(bprimes,ind); bprimes=hellbadprimes(el); return(prod(ind=1,length(bprimes),hellissmooth(el,P,bprimes[ind])))} \\ is a multiple of the index of the everywhere smooth points, probably it is not always the index. \\ if el is minimal this is just the lcm of the local tamagawa-factors hellbindex(el) = { local(ellmin,bprimes,u,ind,produit,pri); produit=1; bprimes=hellbadprimes(el); ellmin=ellchangecurve(el,ellglobalred(el)[2]); u=ellglobalred(el)[2][1]; for(ind=1,length(bprimes), pri=bprimes[ind]; if(u %pri, produit=lcm(elllocalred(el,pri)[4],produit), produit=lcm(pri^(valuation(u,pri)-1) * hellnp(ellmin,pri) *elllocalred(el,pri)[4], produit))); return(produit) } hellsmoothmult(el,P)= { local(mp,d); fordiv(hellbindex(el),d, mp=ellpow(el,P,d); if(hellisallsmooth(el,mp), return([mp,d]);break();)); } \\========================================================= \\ initialisation of an elliptic curve for height calculation \\ hellinit is defined for \\ a 5term vector representing a Weierstrass equation of an elliptic curve over Q and \\ a number bigO (ex. =20,100,1000) being the precision of the formal group law \\ \\ the output is a vector [e,xeris,zeries,sigmacoeff] where \\ e is the elliptic curve \\ xeries is x(t) = t^-2 + \\ zeries is the formal logarithm of the formal group z(t) = t + \\ sigmacoeff is a list of polynomials in Q[a] as they appear in the sigmafunctions associated to alpha = a \\ the "a" here is -3/2 s_2 of Perrin-Riou. But by this it is garanteed \\ that a lies in ZZ_p for p>2. hellinit(el,bigO) = { local(e,s,f,xs,zs,ys,si0,si,sico); e=ellinit(el); \\using formal group law to calculate the x, y and z series f = Mod(t^3,t^bigO) + Mod(e.a1 * t + e.a2 * t^2,t^bigO) * w + Mod(e.a3 + e.a4 * t,t^bigO) * w^2 + Mod(e.a6,t^bigO)* w^3; s = f; while(poldegree(s,w)>0, s=truncate(subst(s,w,f)+O(w)^bigO)); \\ the series ws = (lift(s)+O(t)^bigO); ys = -1 / ws; xs = -t * ys; zs = intformal(1/(2*ys+e.a1*xs+e.a3)*deriv(xs,t),t); \\the sigma functionS si0 = z*exp(a/3*z^2-intformal(intformal(ellwp(e,z)-1/z^2))); \\if we write a/3 here we know that for all p>2, there is a a in ZZ_p si= subst(si0,z,zs); sico=Vec(si); return([e,xs,zs,sico]) } \\hellinitp is defined for \\ a 5term vector representing a Weierstrass equation of an elliptic curve over Q, \\ a prime p and \\ a number bigO (ex. =20,100,1000) being the precision of the formal group law \\ \\ the output is a vector [e,alpha,sigma] where \\ e is the elliptic curve \\ alpha is -3/2 s_2 where s_2 is the second p-adic Eisenstein series. \\ sigma is the canonical p-adic sigma-function, here as a polynomial t+... with coefficients in Z_p, hellinitp(el,p,bigO) = { local(e,s,f,xs,zs,ys,si0,si,sico,al,ind,co,cod,co0,co1); e=ellinit(el); if(hellisgoodordinary(e,p), \\using formal group law to calculate the x, y and z series f = Mod(t^3,t^bigO) + Mod(e.a1 * t + e.a2 * t^2,t^bigO) * w + Mod(e.a3 + e.a4 * t,t^bigO) * w^2 + Mod(e.a6,t^bigO)* w^3; s = f; while(poldegree(s,w)>0, s=truncate(subst(s,w,f)+O(w)^bigO)); \\ the series ws = (lift(s)+O(t)^bigO); ys = -1 / ws; xs = -t * ys; zs = intformal(1/(2*ys+e.a1*xs+e.a3)*deriv(xs,t),t); \\the sigma functionS si0 = z*exp(a/3*z^2-intformal(intformal(ellwp(e,z)-1/z^2))); \\if we write a/3 here we know that for all p>2, there is a a in ZZ_p si= subst(si0,z,zs); sico=Vec(si); \\calculating \alpha al=a; for(ind=3, length(sico), co=subst(sico[ind],a,al); cod=denominator(content(co)); if(cod % p,, co0=Mod(cod*polcoeff(co,0),p); co1=Mod(cod*polcoeff(co,1),p); al=subst(al,a, p* a+lift(-co0/co1)); )); al=subst(al,a,O(p^0)); \\putting it into the sigmafunction \\attention : les coefficients commencent par le t^1 , il faut multiplier par t ! si=t*Polrev(subst(sico,a,al),t); return([e,p,al,si]), print(p" is not an ordinary good prime"); return(0)) } \\hellatp can be used to specify the prime p in a output he of hellinit to get a result like in hellinitp hellatp(he,p) = { local(e,si,sico,al,ind,co,cod,co0,co1); e=he[1]; sico=he[4]; if(hellisgoodordinary(e,p), \\calculating \alpha the "second eisenstein series" al=a; for(ind=3, length(sico), co=subst(sico[ind],a,al); cod=denominator(content(co)); if(cod % p,, co0=Mod(cod*polcoeff(co,0),p); co1=Mod(cod*polcoeff(co,1),p); al=subst(al,a, p* a+lift(-co0/co1)); )); al=subst(al,a,O(p^0)); \\putting it into the sigmafunction \\attention : les coefficients commencent par le t^1 , il faut multiplier par t ! si=t*Polrev(subst(sico,a,al),t); return([e,p,al,si]), print(p" is not an good ordinary prime"); return(0)) } \\========================================================= \\ p-adic height of the point pt. hep is the output of hellinitp or hellatp hellheight(hep,pt)= { local(e,ptp,d,np,p,m,tval,ttt,al,si,sigmaval); \\ check if everything is ok if(hep, e=hep[1]; if(ellisoncurve(e,pt), p=hep[2]; al=hep[3]; si=hep[4]; \\put the point into the subgroup of good reduction everywhere ptp = hellsmoothmult(e,pt); d=ptp[2]; ptp=ptp[1]; \\put the point in the formal group np=hellnp(e,p); while(hellee(ptp) % p != 0, m=divisors(np)[2]; ptp=ellpow(e,ptp,m); d=d*m; np=np/m;); tval=valuation(helltt(ptp),p); \\define precision ttt=helltt(ptp)+O(p^(padicprec(al,p)+2*tval+1)); \\value of the sigma function sigmaval=subst(si,t,ttt); if(sigmaval % p^padicprec(sigmaval,p), return(1/d^2*log(hellee(ptp)/sigmaval)), print("division by zero"); return(0)), print(pt" is not on the elliptic curvre"); return(0)), print(p" is not a ordinary good prime"); return(0))} \\ this should be changed. In fact there is a direct formula not involving the sum of P and Q which should be much faster. hellbil(hep,pt1,pt2) = hellheight(hep,elladd(hep[1],pt1,pt2))-hellheight(hep,pt1)-hellheight(hep,pt2) hellreg(hep,ba) = { local(rank,hmat,hmaten,ind,jnd); rank=length(ba); if(rank<1,return(1), hmat=matrix(rank,rank); for(ind=1,rank, hmat[ind,ind]=2*hellheight(hep,ba[ind]) ); for(ind=1,rank, for(jnd=ind+1,rank, hmaten=hellbil(hep,ba[ind],ba[jnd]); hmat[ind.jnd]=hmaten; hmat[jnd,ind]=hmaten; )); return(matdet(hmat)))} \\========================================================= helleuler(el,ba,p,printoutput=1) = {local(rank,hep,reg,gt,tamanu,lowerbound,lred,lredstr); rank=length(ba); bp=hellbadprimes(el); tamanu=valuation(ellglobalred(el)[3],p); \\fshabound=valuation(round(ellSha(el,ba)),p); gt=valuation(elltors(el)[1],p); lred=elllocalred(el,p); if(lred[2]==1, lredstr=concat(" good",if(hellnp(el,p)%p == 0,", anomalous "," "))); if(lred[2]==1 && hellisgoodordinary(el,p) ==0, lredstr=concat(lredstr,", supersingular ")); if(lred[2]==2, lredstr="II"); if(lred[2]==3, lredstr="III"); if(lred[2]==4, lredstr="IV"); if(lred[2]==-1, lredstr="I_0 *"); if(lred[2]==-2, lredstr="II*"); if(lred[2]==-3, lredstr="III*"); if(lred[2]==-4, lredstr="IV*"); if(lred[2]>4, lredstr=concat("I_",lred[2]-4)); if(lred[2]>4, if(ellap(el,p)==-1, lredstr=concat(lredstr, ", non-"),lredstr=concat(lredstr, ", ")); lredstr=concat(lredstr,"split,")); if(lred[2]<-4, lredstr=concat(concat("I_",- lred[2] -4)," *")); if(lred[2]!=1,lredstr=concat(lredstr,concat(" c = ",lred[4]))); \\output if(printoutput >0, print(" curve of conductor "ellglobalred(el)[1]); print(" reduction at p "lredstr); print(" global torsion "p"^"gt); print(" tamagawa "p"^"tamanu);); \\ordinary case if( hellisgoodordinary(el,p), if(rank>0, hep=hellinitp(el,p,if(p<50,p+10,60)); reg = valuation(hellreg(hep,ba),p), reg= 0 ); lowerbound= reg -rank +tamanu + 2*valuation(hellnp(el,p),p)-2*gt; if(printoutput >0, print(" regulator "p"^"reg); \\print(" analytic Sha(p) "p"^"fshabound); print(""); print(" Euler char/Sha "p"^"lowerbound); ); return([lowerbound,reg]), return([-1,-1]))} \\========================================================= addhelp(hellinit,"Initialisation of the sigma function. Given a 5-term vector el, an integer bigO>4, hellinit(el,bigO) computes a vector [e,xs,zs,sigmacoeff] where e is the elliptic curve el, xs is the formal power series x(t), zs is the formal logarithm z(t) and sigmacorff is a list of coefficients of the sigma function exp(a/3 z^2)*sigma_0 for an unknown parameter a. This function is used for computing p-adic heights with hellatp and hellheight"); addhelp(hellinitp,"Initialisation of the canonical p-adic sigma function. Given a 5-term vector el, a prime p>2 and integer bigO>4, hellinitp(el,p,bigO) computes a vector [e,alpha,sigma] where e is the elliptic curve el and alpha is -3/2 s_2 with s_2 is the second p-adic Eisenstein series. Finally sigma is the canonical p-adic sigma-function, here as a polynomial t+... with coefficients in Z_p. This function is used for computing p-adic heights with hellheight"); addhelp(hellatp,"Initialisation of the canonical p-adic sigma function. Given the output he of the function hellinit and a prime p>2, hellinitp(el,p,bigO) computes a vector [e,alpha,sigma] where e is the elliptic curve el and alpha is -3/2 s_2 with s_2 is the second p-adic Eisenstein series. Finally sigma is the canonical p-adic sigma-function, here as a polynomial t+... with coefficients in Z_p. This function is used for computing p-adic heights with hellheight"); addhelp(hellheight,"Given the output hep of either hellinitp or hellatp and a point P, the function hellheight(hep,P) computes the canonical p-adic height of the point P"); addhelp(hellbil,"Given the output hep of either hellinitp or hellatp and two points P and Q, the function hellheight(hep,P,Q) computes the canonical p-adic height pairing of the point P"); addhelp(hellreg,"Given the output hep of either hellinitp or hellatp and a basis ba of the Mordell-Weil modulo torsion, the function hellreg(hep,ba) computes the canonical p-adic regulator"); addhelp(hellisgoodordinary,"hellisgoodordinary(el,p) checks if the elliptic curve el has good and ordinary reduction at the prime p"); addhelp(hellisanomalous,"hellisanomalous(el,p) checks if the elliptic curve el has anomalous reduction at the prime p "); addhelp(hellnp,"hellnp(el,p) computes the number of non-singular points on the reduction of the elliptic curve el at the prime p"); addhelp(hellissmooth,"hellissmooth(el,P,p) checks if the point P on the elliptic curve e reduces to a non-singular point in the reduction at the prime p"); addhelp(hellisallsmooth,"hellallissmooth(el,P) checks if the point P on the elliptic curve e reduces to a non-singular point at all primes, i.e. if it belongs to the connected component of the Néron model"); addhelp(helleuler,"helleuler(el,ba,p) computes several invariants appearing in the euler-characteristic of the elliptic curve el at the prime p. ba must be a basis of the Mordell-Weil group modulo torsion. If the reduction is good and ordinary, the result is a vector containing the valuation of the euler-char (up to Sha) and the valuation of the regulator.");