| Download
All published worksheets from http://sagenb.org
Project: sagenb.org published worksheets
Views: 168703Image: ubuntu2004
from sage.all import * class AntiCyclotomicTower: def __init__(self,E,D,p,prec=2000): """ E : elliptic curve D : quadratic imaginary number field satisfying Heegner hypothesis p : good ordinary prime prec : bits of precision for complex-valued calculations """ self._E = E self._D = D self._p = p self._prec = prec def tau(self, f): p = self._p sqrtD = self._sqrtD c = self._c d = c*sqrtD A,B,_ = f return (-B + d)/(2*A) def qf_to_ideal(self, f): c = self._c OK = self._OK sqrtD = self._sqrtD (A,B,C) = f if A%c == 0: A,C=C,A return OK.ideal([A,(-B+c*sqrtD)/2]) def mult(self, alpha, beta): v_alpha = self._v_alpha v_P1 = self._v_P1 return v_alpha[v_P1.index(to_P1(alpha*beta))] def mult_order_bq(self,alpha): n = 1 id = self._e.reduced_form() x = alpha.reduced_form() while x!= id: x = (x*alpha).reduced_form() n +=1 return n def z_point_bq(self,F,list_of_qf): self._w_Q = prime_to_p_bq = self._prime_to_p_bq w_tau = [self.tau(f) for f in list_of_qf] C = F.base_field() w_tau_numerical = [C(t) for t in w_tau] phi = self._E.modular_parametrization() points = [F.point(phi(t),check=False) for t in w_tau_numerical] return sum(points) def padic_height(self,n=2, check=False): """ Let K = Q(sqrt(D)) and c = p^2. The nth layer ``M_{n-1}'' of the anticyclotomic tower sits inside K_{p^n} | M_{n-1} | K | Q Let n = 2. Then [K_{p^2}:M_1] = p pm 1, with + 1 if p splits in K and -1 if p is inert in K. This returns the p-adic height of a point P in M_{n-1} EXAMPLES: sage: E = EllipticCurve('57a') sage: D = -8 sage: p = 5 sage: prec = 2000 sage: T = AntiCyclotomicTower(E,D,p,prec) sage: T.padic_height(2) 2 + 2*5 + 2*5^2 + 5^4 + 4*5^5 + 4*5^6 + 3*5^7 + 4*5^8 + 5^9 + O(5^12) sage: E = EllipticCurve('57a') sage: p = 5 sage: D = -56 sage: prec = 2000 sage: T = AntiCyclotomicTower(E,D,p,prec) sage: T.padic_height(2) 3 + 2*5 + 5^2 + 4*5^5 + 2*5^6 + 4*5^8 + 5^9 + 4*5^10 + 4*5^11 + O(5^12) sage: E = EllipticCurve('89a') sage: p = 5 sage: prec = 2000 sage: D = -8 sage: T = AntiCyclotomicTower(E,D,p,prec) sage: T.padic_height(2) 4 + 3*5 + 4*5^2 + 5^3 + 4*5^4 + 5^5 + 5^7 + 4*5^8 + 4*5^9 + 5^10 + 4*5^11 + O(5^12) sage: D = -67 sage: T = AntiCyclotomicTower(E,D,p,prec) sage: T.padic_height(2) 3 + 2*5 + 5^2 + 5^3 + 2*5^4 + 2*5^6 + 4*5^7 + 5^8 + 3*5^9 + 2*5^10 + 2*5^11 + O(5^12) sage: E = EllipticCurve('158b1') sage: D = -7 sage: p = 5 sage: prec = 2000 sage: T = AntiCyclotomicTower(E,D,p,prec) sage: T.padic_height(2) 3*5 + 4*5^2 + 2*5^3 + 4*5^4 + 4*5^5 + 3*5^6 + 3*5^7 + 2*5^8 + 2*5^9 + 2*5^10 + O(5^13) sage: E = EllipticCurve('214b1') sage: D = -7 sage: p = 5 sage: prec = 2000 sage: T = AntiCyclotomicTower(E,D,p,prec) sage: T.padic_height(2) 2 + 2*5 + 3*5^2 + 2*5^3 + 2*5^4 + 2*5^5 + 3*5^6 + 2*5^7 + 2*5^8 + 2*5^9 + 4*5^10 + 2*5^11 + O(5^13) sage: E = EllipticCurve('158b1') sage: D = -7 sage: p = 5 sage: prec = 6000 sage: T = AntiCyclotomicTower(E,D,p,prec) sage: time T.padic_height(3) 3*5 + 2*5^2 + 3*5^3 + 2*5^4 + 5^5 + 4*5^6 + 3*5^7 + 4*5^8 + O(5^9) """ proof.number_field(False) prec = self._prec if n <= 2: assert prec >= 2000 elif n == 3: assert prec >= 5000 else: print "Warning: make sure prec is large enough!" E = self._E D = self._D p = self._p self._n = n self._c =c = p**n K = QuadraticField(D, names='sqrtD') if len(factor(p*K)) == 2: print "Warning: p splits!" self._sqrtD = sqrtD = K.gen() OK = self._OK = K.ring_of_integers() H = heegner_points(E.conductor(), D, c) beta = H.betas()[0] self._QF = QF = [x.quadratic_form() for x in H.points(beta)] QF1 = [a for a in QF if len(self.qf_to_ideal(a).gens_reduced()) == 1] for i in QF1: if (i*i).reduced_form() == i.reduced_form() and (i*i*i).reduced_form() == i.reduced_form(): self._e = e = i break self._bq_ord = bq_ord = [self.mult_order_bq(x) for x in QF] prime_to_p_bq = [] ord_p_bq = [] for i in range(len(QF)): x = bq_ord[i] if x%p != 0: prime_to_p_bq.append(QF[i]) if x == p: ord_p_bq.append(QF[i]) self._prime_to_p_bq = prime_to_p_bq self._ord_p_bq = ord_p_bq ord5 = ord_p_bq[0] coset1 = [] coset2 = [] coset3 = [] coset4 = [] coset5 = [] self._bq_red = bq_red = [x.reduced_form() for x in QF] for g in prime_to_p_bq: coset1.append(QF[bq_red.index( (ord5*g).reduced_form() )]) coset2.append(QF[bq_red.index((ord5*ord5*g).reduced_form())]) coset3.append(QF[bq_red.index((ord5*ord5*ord5*g).reduced_form())]) coset4.append(QF[bq_red.index((ord5*ord5*ord5*ord5*g).reduced_form())]) coset5.append(QF[bq_red.index((ord5*ord5*ord5*ord5*ord5*g).reduced_form())]) self._w = w = self.z_point_bq(self._E.change_ring(ComplexField(self._prec)), prime_to_p_bq) print "Here are the conjugates of the Heegner points" print "x(z_1) = ", CC(w[0]) print "y(z_1) = ", CC(w[1]) sigmaw = self.z_point_bq(self._E.change_ring(ComplexField(self._prec//10)), coset1) print "x(sigma(z_1)) = ", CC(sigmaw[0]) print "y(sigma(z_1)) = ", CC(sigmaw[1]) sigma2w = self.z_point_bq(self._E.change_ring(ComplexField(self._prec//10)), coset2) print "x(sigma^2(z_1))= ", CC(sigma2w[0]) print "y(sigma^2(z_1))= ", CC(sigma2w[1]) sigma3w = self.z_point_bq(self._E.change_ring(ComplexField(self._prec//10)), coset3) print "x(sigma^3(z_1)) = ", CC(sigma3w[0]) print "y(sigma^3(z_1)) = ", CC(sigma3w[1]) sigma4w = self.z_point_bq(self._E.change_ring(ComplexField(self._prec//10)), coset4) print "x(sigma^4(z_1)) = ", CC(sigma4w[0]) print "y(sigma^4(z_1)) = ", CC(sigma4w[1]) sigma5w = self.z_point_bq(self._E.change_ring(ComplexField(self._prec//10)), coset5) print "x(sigma^5(z_1)) = ", CC(sigma5w[0]) print "y(sigma^5(z_1)) = ", CC(sigma5w[1]) #Enew = E.change_ring(ComplexField(53)) #print "trace = ", Enew(sigmaw) + Enew(sigma2w) + Enew(sigma3w) + Enew(sigma4w) + Enew(w) f = algebraic_dependency(w[0],p**(n-1)) C = ComplexField(self._prec)['x'] all_roots = [i[0] for i in (C(f)).roots()] real_root = [] complex_roots = [] for j in all_roots: try: real_root = [RealField(self._prec)(j)] except TypeError: if j.conjugate() not in complex_roots: complex_roots.append(j) self._complex_roots = complex_roots degf = f.degree() lead_coeff = f[degf] digits = log(RR(lead_coeff))/log(RR(10)) print "------------------------------" print "minpoly of x-coordinates of Heegner points = ", f roots = (GF(self._p)['x'])(f).roots() #print "roots = ", roots r = roots[0][0] ord = E.change_ring(GF(self._p)).lift_x(r).additive_order() ord = ord*prod(E.tamagawa_numbers()) #print "ord = ", ord f = f.monic() L = NumberField(f,names='b') b = L.gen() if floor(digits)//2 < 20: padic_prec = 20 else: padic_prec = floor(digits)//2 F = Qp(self._p,padic_prec) T = QQ['x'] x = T.gen() for i in range(self._p): try: G = F.extension(T(f(x+i)),names='alpha') offset = i if offset == i: break except NotImplementedError: pass G = F.extension(T(f(x+offset)),names='alpha') alpha = G.gen() EG = E.change_ring(G) R = G['x','y'] x,y = R.gens() g = R(E.defining_polynomial()(x,y,1)) S = G['y'] y = S.gen() h = S(g(alpha+offset,y)) hprime = S(h.derivative(y)) sq = h[1]**2 - 4*h[0]*h[2] dd = (F(sq[0])).sqrt() for i in range(20): dd = (dd + sq/dd)/ZZ(2) ycoord = G((-h[1]+dd)/(2*h[2])) P = EG(alpha+offset,ycoord) Q = P for i in range(ord-1): Q = P + Q val_for_sig = -Q[0]/Q[1] sig = E.padic_sigma(self._p,self._p*padic_prec) H = G[['t']] sig = H(sig) sval = sig(val_for_sig) if ord%2 == 0: fn = E.division_polynomial(ord, two_torsion_multiplicity =1) fnval = fn(P[0],P[1]) else: fn = E.division_polynomial(ord) fnval = fn(alpha+offset) h1 = ZZ(1)/ZZ(self._p*ord**2)*(log( (sval.norm())**2/(fnval.norm())**2 )-ord**2*log(F(lead_coeff))) print "height of z_1:" return h1 def height_point(self,pt): n = self._n E = self._E prec = self._prec EC = E.change_ring(ComplexField(prec)) f = algebraic_dependency(pt[0],self._p**(n-1)) print f.discriminant().valuation(5) print f.discriminant().valuation(7) degf = f.degree() lead_coeff = f[degf] digits = log(RR(lead_coeff))/log(RR(10)) roots = (GF(self._p)['x'])(f).roots() print "roots = ", roots r = roots[0][0] ord = E.change_ring(GF(self._p)).lift_x(r).additive_order() ord = ord*prod(E.tamagawa_numbers()) f = f.monic() if floor(digits)//2 < 20: padic_prec = 20 else: padic_prec = floor(digits)//2 F = Qp(self._p,padic_prec) T = QQ['x'] x = T.gen() for i in range(self._p): try: G = F.extension(T(f(x+i)),names='alpha') offset = i if offset == i: break except NotImplementedError: pass G = F.extension(T(f(x+offset)),names='alpha') alpha = G.gen() EG = E.change_ring(G) R = G['x','y'] x,y = R.gens() g = R(E.defining_polynomial()(x,y,1)) S = G['y'] y = S.gen() h = S(g(alpha+offset,y)) hprime = S(h.derivative(y)) sq = h[1]**2 - 4*h[0]*h[2] dd = (F(sq[0])).sqrt() for i in range(20): dd = (dd + sq/dd)/ZZ(2) ycoord = G((-h[1]+dd)/(2*h[2])) P = EG(alpha+offset,ycoord) Q = P for i in range(ord-1): Q = P + Q val_for_sig = -Q[0]/Q[1] sig = E.padic_sigma(self._p,self._p*padic_prec) H = G[['t']] sig = H(sig) sval = sig(val_for_sig) if ord%2 == 0: fn = E.division_polynomial(ord, two_torsion_multiplicity =1) fnval = fn(P[0],P[1]) else: fn = E.division_polynomial(ord) fnval = fn(alpha+offset) h1 = ZZ(1)/ZZ(self._p*ord**2)*(log( (sval.norm())**2/(fnval.norm())**2 )-ord**2*log(F(lead_coeff))) print "h1 = ", h1 ######################### ord = 2*ord Q = P for i in range(ord-1): Q = P + Q val_for_sig = -Q[0]/Q[1] sval = sig(val_for_sig) fn = E.division_polynomial(ord, two_torsion_multiplicity =1) fnval = fn(P[0],P[1]) h2 = ZZ(1)/ZZ(self._p*ord**2)*(log( (sval.norm())**2/(fnval.norm())**2 )-ord**2*log(F(lead_coeff))) print "h2 = ", h2 return h1 def height_paired_against_conjugates(self): n = self._n complex_roots = self._complex_roots E = self._E prec = self._prec EC = E.change_ring(ComplexField(prec)) heights = [] for root in complex_roots: print "for %s and its conjugate "%CC(root) pt = EC.lift_x(root) + EC.lift_x(root.conjugate()) try: RR(pt[0]), RR(pt[1]) except TypeError: print "ooops..." f = algebraic_dependency(pt[0],p**(n-1)) degf = f.degree() lead_coeff = f[degf] digits = log(RR(lead_coeff))/log(RR(10)) roots = (GF(self._p)['x'])(f).roots() #print "roots = ", roots r = roots[0][0] ord = E.change_ring(GF(self._p)).lift_x(r).additive_order() ord = ord*prod(E.tamagawa_numbers()) #print "ord = ", ord f = f.monic() if floor(digits)//2 < 20: padic_prec = 20 else: padic_prec = floor(digits)//2 F = Qp(self._p,padic_prec) T = QQ['x'] x = T.gen() for i in range(self._p): try: G = F.extension(T(f(x+i)),names='alpha') offset = i if offset == i: break except NotImplementedError: pass G = F.extension(T(f(x+offset)),names='alpha') alpha = G.gen() EG = E.change_ring(G) R = G['x','y'] x,y = R.gens() g = R(E.defining_polynomial()(x,y,1)) S = G['y'] y = S.gen() h = S(g(alpha+offset,y)) hprime = S(h.derivative(y)) sq = h[1]**2 - 4*h[0]*h[2] dd = (F(sq[0])).sqrt() for i in range(20): dd = (dd + sq/dd)/ZZ(2) ycoord = G((-h[1]+dd)/(2*h[2])) P = EG(alpha+offset,ycoord) Q = P for i in range(ord-1): Q = P + Q val_for_sig = -Q[0]/Q[1] sig = E.padic_sigma(self._p,self._p*padic_prec) H = G[['t']] sig = H(sig) sval = sig(val_for_sig) if ord%2 == 0: fn = E.division_polynomial(ord, two_torsion_multiplicity =1) fnval = fn(P[0],P[1]) else: fn = E.division_polynomial(ord) fnval = fn(alpha+offset) h1 = ZZ(1)/ZZ(self._p*ord**2)*(log( (sval.norm())**2/(fnval.norm())**2 )-ord**2*log(F(lead_coeff))) print h1 heights.append(h1) print "----------------------" return heights def z_points(self, alpha, F): alpha_p = self._alpha_p v = [] for i in range(1): v.append(sum(self.z_point(alpha, F))) #print "v = ", v alpha = self.mult(alpha, alpha_p) return v
E = EllipticCurve('57a') D = -8 p = 5 prec = 2000 T = AntiCyclotomicTower(E,D,p,prec) T.padic_height(2)
Here are the conjugates of the Heegner points
x(z_1) = 1.09134357351891 - 9.63652225790997e-596*I
y(z_1) = -0.919649689611060 - 1.09577255026657e-595*I
x(sigma(z_1)) = 1.28240225474401 - 0.182500350994469*I
y(sigma(z_1)) = -0.761690770112933 + 0.117006496908598*I
x(sigma^2(z_1))= 1.67723875767367 + 0.0866463691344989*I
y(sigma^2(z_1))= -1.39041234698688 - 0.149731706982934*I
x(sigma^3(z_1)) = 1.67723875767367 - 0.0866463691344989*I
y(sigma^3(z_1)) = -1.39041234698688 + 0.149731706982934*I
x(sigma^4(z_1)) = 1.28240225474401 + 0.182500350994469*I
y(sigma^4(z_1)) = -0.761690770112933 - 0.117006496908598*I
x(sigma^5(z_1)) = 1.09134357351891 - 5.27392019267801e-54*I
y(sigma^5(z_1)) = -0.919649689611060 - 9.71087656786162e-54*I
------------------------------
minpoly of x-coordinates of Heegner points = 18034072681*x^5 - 126430131580*x^4 + 352783410220*x^3 - 489834319200*x^2 + 338504989540*x - 93144838864
height of z_1:
2 + 2*5 + 2*5^2 + 5^4 + 4*5^5 + 4*5^6 + 3*5^7 + 4*5^8 + 5^9 + 5^12 + 5^13 + 3*5^14 + 2*5^15 + 4*5^16 + 5^17 + O(5^18)
R.<x> = CC['x'] f = 18034072681*x^5 - 126430131580*x^4 + 352783410220*x^3 -489834319200*x^2 + 338504989540*x - 93144838864 roots = [x[0] for x in f.roots()]; roots
[1.09134357351899, 1.28240225474406 - 0.182500350994335*I, 1.28240225474406 + 0.182500350994335*I, 1.67723875767358 - 0.0866463691342926*I, 1.67723875767358 + 0.0866463691342926*I]
Here we check that we're choosing the right y-coordinates of the conjugates
-E.lift_x(roots[0])
(1.09134357351899 : -0.919649689611003 : 1)
-E.lift_x(roots[1])
(1.28240225474406 - 0.182500350994335*I : -0.761690770113188 + 0.117006496908515*I : 1)
-E.lift_x(roots[2])
(1.28240225474406 + 0.182500350994335*I : -0.761690770113188 - 0.117006496908515*I : 1)
-E.lift_x(roots[3])
(1.67723875767358 - 0.0866463691342926*I : -1.39041234698675 + 0.149731706982566*I : 1)
-E.lift_x(roots[4])
(1.67723875767358 + 0.0866463691342926*I : -1.39041234698675 - 0.149731706982566*I : 1)
So this is what the trace of z1 computed directly by hand looks like
sum(-E.lift_x(i) for i in roots) #this is tr z_1, which is 3*z0
(3.99999999999904 - 1.37667655053519e-14*I : -6.99999999999718 + 4.13002965160558e-14*I : 1)
Bertolini's formula says that we should get trz1=(ap−ap−1∗(1+p))(z0)=−1∗z0
ap = E.ap(5); ap
-3
ap - ap^-1*(1+p)
-1
z0 = E.heegner_point(D).point_exact(); z0
(1 : 0 : 1)
-z0
(1 : -1 : 1)
3*z0
(4 : -7 : 1)