2011-07-01-saturate

327 days ago by arabindr

K.<a> = NumberField(x^2 - x - 1) import psage.ellcurve.lseries.sqrt5 as sqrt5w 
       
class ReductionMap: """ INPUT: - E -- Elliptic curve over K=Q(sqrt(5)) - P -- a prime in Q that splits in K - ell -- a prime that exactly divides the cardinality of E mod P. OUTPUT: - a callable map from E(K) to GF(ell) """ def __init__(self, E, P, ell): self.E = E self.P = P self.ell = ell self.v = K.factor(self.P) Emod0 = E.change_ring(self.v[0][0].residue_field()) Emod1 = E.change_ring(self.v[1][0].residue_field()) self.Emod0 = Emod0 self.Emod1 = Emod1 self.GFell = GF(ell) # compute the ell-torsion subgroup T0 = Emod0(0).division_points(ell) T1 = Emod1(0).division_points(ell) if len(T0) == ell: self.T = T0 self.Emod = Emod0 elif len(T1) == ell: self.T = T1 self.Emod = Emod1 else: raise ValueError, "must have Emod[ell] of order ell" # compute multiplier self.card = self.Emod.cardinality() self.multiplier = ZZ(self.card / ell) # make discrete log table # 1. find first nonzero element g in T i = 0 g = self.T[i] while g == 0: i += 1 g = self.T[i] dlog = {g:1} # 2. set h to 2*g, 3*g, .., recording the multiple h = g + g for n in range(2,ell+1): dlog[h] = n h += g self.dlog = dlog def __call__(self, Q): # Reduce Q to the curve over the finite field. # Then multiple the result by #Emod/ell. # Finally write in terms of generator using dlog table. Qbar = self.Emod(Q) n = self.dlog[self.multiplier * Qbar] return self.GFell(n) 
       
def nice_prime_finder(E, norm_bound, ell): prime_set = prime_range(norm_bound) v=[] t = 0 for l in prime_set: if l!=5 and len(K.factor(l))==2: for P in K.primes_above(l): n = E.change_ring(P.residue_field()).cardinality() if n%ell ==0 and (n/ell)%ell != 0: v.append(l) return v 
       
def saturate_prime(E, Q, ell): norm_bound = 100 w=[] while len(w) < 5: w+=nice_prime_finder(E, norm_bound, ell) norm_bound += 100 psi0 = ReductionMap(E, w[0], ell) psi1 = ReductionMap(E, w[1], ell) i = 2 M = matrix(GF(ell), [[psi0(Q[0]), psi0(Q[1])], [psi1(Q[0]),psi1(Q[1])]]) while M[0] == M[1]: psii = ReductionMap(E, w[i], ell) M = matrix(GF(ell), [[psi0(Q[0]), psi0(Q[1])], [psii(Q[0]),psii(Q[1])]]) i+=1 while det(M) == 0: a1 = ZZ(M[0][0]) a2 = ZZ(M[0][1]) if a1 == a2 == 0: a1 = ZZ(M[1][0]) a2 = ZZ(M[1][1]) g = gcd(a1,a2) a1/=g a2/=g Q2 = (-ZZ(a2)*E(Q[0]) + ZZ(a1)*E(Q[1])).division_points(ell)[0] if a1 == 0: Q[0] = Q2 else: Q[1] = Q2 i = 2 M = matrix(GF(ell), [[psi0(Q[0]), psi0(Q[1])], [psi1(Q[0]),psi1(Q[1])]]) while M[0] == M[1]: psii = ReductionMap(E, w[i], ell) M = matrix(GF(ell), [[psi0(Q[0]), psi0(Q[1])], [psii(Q[0]),psii(Q[1])]]) i+=1 return Q 
       
E = EllipticCurve([0,-a,1,-a-1,2*a+1] ) Q = E.gens(); print 'Q is', Q; print Q[0].height(), Q[1].height() 
       
Q is [(0 : -a - 1 : 1), (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)]
0.166157999281333 1.94042291278625
Q is [(0 : -a - 1 : 1), (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)]
0.166157999281333 1.94042291278625
print saturate_prime(E, Q, 3) print Q[0].height(),Q[1].height() 
       
[(0 : -a - 1 : 1), (a + 1 : -a - 1 : 1)]
0.166157999281333 0.215602545865140
[(0 : -a - 1 : 1), (a + 1 : -a - 1 : 1)]
0.166157999281333 0.215602545865140
def real_periods(E): embs = K.embeddings(RR) v=[1,1] i = 0 while i < 2: t = E.period_lattice(embs[i]).basis() if real(t[1])==0: v[i] = 2*t[0] else: v[i] = t[0] i+=1 return v[0]*v[1] 
       
def hp(X, Y): return 0.5((X+Y).height() - X.height() - Y.height()) 
       
reg = hp(Q[0], Q[0])*hp(Q[1],Q[1]) - hp(Q[1],Q[0])^2 
       
Traceback (click to the left of this block for traceback)
...
TypeError: 'sage.rings.real_mpfr.RealLiteral' object is not callable
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_88.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("cmVnID0gaHAoUVswXSwgUVswXSkqaHAoUVsxXSxRWzFdKSAtIGhwKFFbMV0sUVswXSleMg=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpsl0goS/___code___.py", line 3, in <module>
    exec compile(u'reg = hp(Q[_sage_const_0 ], Q[_sage_const_0 ])*hp(Q[_sage_const_1 ],Q[_sage_const_1 ]) - hp(Q[_sage_const_1 ],Q[_sage_const_0 ])**_sage_const_2 
  File "", line 1, in <module>
    
  File "/tmp/tmpP4LMdT/___code___.py", line 4, in hp
    return _sage_const_0p5 ((X+Y).height() - X.height() - Y.height())
TypeError: 'sage.rings.real_mpfr.RealLiteral' object is not callable
def conjectural_sha(E, reg, Lstar=1): M = E.tamagawa_product_bsd() sha = RR(sqrt(5))*Lstar*(E.torsion_order())^2/((real_periods(E))*reg*M) return sha 
       
print conjectural_sha(E, reg, 2.88288222151816) # This is confusing me. I thought the problem was with the regulator of points, but my rewrite changes nothing. Help? 
       
7.99999999999992
7.99999999999992