Copyofrank_two_computation_sha_reg_MW_for_Ashwath

323 days ago by Ben.LeVeque

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_lower, norm_bound, ell): prime_set = prime_range(norm_lower, 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_old(E, Q, ell): norm_bound = 100 w=[] while len(w) < 5: w+=nice_prime_finder(E, norm_bound - 100, 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 try: Q2 = (ZZ(a2)*E(Q[0]) - ZZ(a1)*E(Q[1])).division_points(ell)[0] except IndexError: break if a1: 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 
       
def saturate_prime(E, Q, ell): w=[] ind = False norm_bound = 100 while ind == False: c = True t = len(w) while len(w) < 15: w+=nice_prime_finder(E, norm_bound - 100, 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 print w while det(M) == 0: print ell print M 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 print a1,a2 print Q[0], Q[1] try: Q2 = (ZZ(a2)*E(Q[0]) - ZZ(a1)*E(Q[1])).division_points(ell)[0] except IndexError: break if a1 == 0: Q[1] = Q2 else: Q[0] = Q2 M = matrix(GF(ell), [[psi0(Q[0]), psi0(Q[1])], [psi1(Q[0]),psi1(Q[1])]]) if det(M): ind = True 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
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 conjectural_sha(E, Lstar): def hp(X, Y): return (X+Y).height() - X.height() - Y.height() Q = E.gens() for ell in prime_range(7): Q = saturate_prime(E, Q, ell) reg = (hp(Q[0], Q[0])*hp(Q[1],Q[1]) - hp(Q[1],Q[0])^2) M = E.tamagawa_product_bsd() sha = RR(sqrt(5))*Lstar*(E.torsion_order())^2/((real_periods(E))*reg*M) return sha, reg, Q 
       
shoutput = conjectural_sha(E, 1.44144111075908) 
       
WARNING: Output truncated!  
full_output.txt



[19, 59, 89, 139, 151, 229, 269, 269, 281, 311, 379, 389, 401, 419, 421,
421, 431, 439, 439, 491, 499, 499]
[19, 29, 59, 61, 131, 131, 139, 151, 179, 191, 239, 281, 349, 359, 401,
401, 409, 431, 449, 449, 461, 479, 499]
3
[2 0]
[0 0]
1 0
(0 : -a - 1 : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]

...

3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
Traceback (click to the left of this block for traceback)
...
__SAGE__
WARNING: Output truncated!  
full_output.txt



[19, 59, 89, 139, 151, 229, 269, 269, 281, 311, 379, 389, 401, 419, 421, 421, 431, 439, 439, 491, 499, 499]
[19, 29, 59, 61, 131, 131, 139, 151, 179, 191, 239, 281, 349, 359, 401, 401, 409, 431, 449, 449, 461, 479, 499]
3
[2 0]
[0 0]
1 0
(0 : -a - 1 : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]

...

3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
3
[0 0]
[2 0]
1 0
(a + 1 : a : 1) (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_210.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("c2hvdXRwdXQgPSBjb25qZWN0dXJhbF9zaGEoRSwgMS40NDE0NDExMTA3NTkwOCk="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpj85YKG/___code___.py", line 3, in <module>
    exec compile(u'shoutput = conjectural_sha(E, _sage_const_1p44144111075908 )
  File "", line 1, in <module>
    
  File "/tmp/tmphUQqfj/___code___.py", line 8, in conjectural_sha
    Q = saturate_prime(E, Q, ell)    
  File "/tmp/tmptAsA8H/___code___.py", line 36, in saturate_prime
    Q2 = (ZZ(a2)*E(Q[_sage_const_0 ]) - ZZ(a1)*E(Q[_sage_const_1 ])).division_points(ell)[_sage_const_0 ]
  File "/home/bleveque/sage/local/lib/python2.6/site-packages/sage/schemes/elliptic_curves/ell_point.py", line 1005, in division_points
    if P == 0:
  File "element.pyx", line 886, in sage.structure.element.Element.__richcmp__ (sage/structure/element.c:7306)
  File "element.pyx", line 831, in sage.structure.element.Element._richcmp (sage/structure/element.c:6757)
  File "element.pyx", line 813, in sage.structure.element.Element._richcmp_ (sage/structure/element.c:6502)
  File "element.pyx", line 858, in sage.structure.element.Element._richcmp (sage/structure/element.c:7183)
  File "/home/bleveque/sage/local/lib/python2.6/site-packages/sage/schemes/elliptic_curves/ell_point.py", line 373, in __cmp__
    def __cmp__(self, other):
KeyboardInterrupt

__SAGE__
print shoutput 
       
(1.00000000000001, 0.0853096123084668, [(0 : -a - 1 : 1), (a + 1 : -a -
1 : 1)])
(1.00000000000001, 0.0853096123084668, [(0 : -a - 1 : 1), (a + 1 : -a - 1 : 1)])
f = open('/home/bleveque/updated_list_copy.txt') discard = f.readline() g = open('/home/bleveque/conjectural_sha_results.txt','w') 
       
def ainv_grabber(line): return eval('['+line.split('[')[1].split(']')[0]+']') 
       
def rank_grabber(line): return str('('+line.split('(')[1].split(')')[0]+')') 
       
f.readline().split()[-4] 
       
'0.359928959498'
'0.359928959498'
def lstar_grabber(line): return eval(line.split()[-4]) 
       
rank_grabber(f.readline()) 
       
'(0, 0)'
'(0, 0)'
sha_list = [] reg_list = [] Q_list = [] for line in f.readlines(): if rank_grabber(line) == '(2, 2)': try: F = EllipticCurve(K,ainv_grabber(line)) print F.a_invariants() shaoutput = list(conjectural_sha(F,lstar_grabber(line))) sha_list.append(shaoutput[0]) reg_list.append(shaoutput[1]) Q_list.append(shaoutput[2]) print shaoutput except ArithmeticError: pass 
       
(0, a - 1, 1, a - 2, -2*a + 3)
3
[1.00000000000063, 0.0853096123084681, [(0 : -a + 1 : 1), (-a + 2 : -a +
1 : 1)]]
(1, -a + 1, a, -a, 0)
[0.999999999999290, 0.105081885636517, [(3/4*a - 3/4 : -a + 9/8 : 1), (0
: -a : 1)]]
(a, 0, 0, a - 3, -a + 2)
3
[1.00000000000118, 0.0929982882741012, [(-1 : 2 : 1), (-a + 1 : a : 1)]]
(a + 1, -a - 1, a, -2*a - 2, a + 1)
3
[1.00000000000262, 0.0989318653886299, [(1 : -a : 1), (0 : -a - 1 : 1)]]
(a, 1, a, -1, 0)
[1.00000000000194, 0.118300192875220, [(-1 : 1 : 1), (0 : -a : 1)]]
(a, 0, 1, -3*a - 2, 2*a + 1)
(a, a + 1, 0, a, 0)
(a, -a - 1, 0, -1, a + 1)
[1.00000000000005, 0.113330086069452, [(0 : -a : 1), (1 : -a : 1)]]
(a + 1, -a, 1, -2, 0)
3
[1.00000000000097, 0.117644186546780, [(-1 : 1 : 1), (1 : -1 : 1)]]
(a + 1, a - 1, a + 1, -a, 0)
(1, -1, 0, a - 2, -a + 2)
3
[1.00000000000237, 0.118355297678359, [(a - 1 : -a + 1 : 1), (-a + 1 : a
- 1 : 1)]]
(a + 1, 0, 1, a - 2, 0)
3
Traceback (click to the left of this block for traceback)
...
IndexError: list index out of range
(0, a - 1, 1, a - 2, -2*a + 3)
3
[1.00000000000063, 0.0853096123084681, [(0 : -a + 1 : 1), (-a + 2 : -a + 1 : 1)]]
(1, -a + 1, a, -a, 0)
[0.999999999999290, 0.105081885636517, [(3/4*a - 3/4 : -a + 9/8 : 1), (0 : -a : 1)]]
(a, 0, 0, a - 3, -a + 2)
3
[1.00000000000118, 0.0929982882741012, [(-1 : 2 : 1), (-a + 1 : a : 1)]]
(a + 1, -a - 1, a, -2*a - 2, a + 1)
3
[1.00000000000262, 0.0989318653886299, [(1 : -a : 1), (0 : -a - 1 : 1)]]
(a, 1, a, -1, 0)
[1.00000000000194, 0.118300192875220, [(-1 : 1 : 1), (0 : -a : 1)]]
(a, 0, 1, -3*a - 2, 2*a + 1)
(a, a + 1, 0, a, 0)
(a, -a - 1, 0, -1, a + 1)
[1.00000000000005, 0.113330086069452, [(0 : -a : 1), (1 : -a : 1)]]
(a + 1, -a, 1, -2, 0)
3
[1.00000000000097, 0.117644186546780, [(-1 : 1 : 1), (1 : -1 : 1)]]
(a + 1, a - 1, a + 1, -a, 0)
(1, -1, 0, a - 2, -a + 2)
3
[1.00000000000237, 0.118355297678359, [(a - 1 : -a + 1 : 1), (-a + 1 : a - 1 : 1)]]
(a + 1, 0, 1, a - 2, 0)
3
Traceback (most recent call last):            try:
  File "", line 1, in <module>
    
  File "/tmp/tmprdCgW_/___code___.py", line 7, in <module>
    exec compile(u"for line in f.readlines():\n    if rank_grabber(line) == '(2, 2)':\n        try:\n            F = EllipticCurve(K,ainv_grabber(line))\n            print F.a_invariants()\n            shaoutput = list(conjectural_sha(F,lstar_grabber(line)))        \n            sha_list.append(shaoutput[_sage_const_0 ])\n            reg_list.append(shaoutput[_sage_const_1 ])\n            Q_list.append(shaoutput[_sage_const_2 ])\n            print shaoutput\n        except ArithmeticError:\n            pass" + '\n', '', 'single')
  File "", line 6, in <module>
    
  File "/tmp/tmpmq3zeu/___code___.py", line 8, in conjectural_sha
    Q = saturate_prime(E, Q, ell)    
  File "/tmp/tmp_y0aBu/___code___.py", line 27, in saturate_prime
    Q2 = (-ZZ(a2)*E(Q[_sage_const_0 ]) + ZZ(a1)*E(Q[_sage_const_1 ])).division_points(ell)[_sage_const_0 ]
IndexError: list index out of range
conjectural_sha(EllipticCurve(K, [a + 1, 0, 1, a - 2, 0]), 1) 
       
2
3
[2 1]
[1 2]
2 1
(0 : 0 : 1) (-1 : 1 : 1)
Traceback (click to the left of this block for traceback)
...
IndexError: list index out of range
2
3
[2 1]
[1 2]
2 1
(0 : 0 : 1) (-1 : 1 : 1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_101.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("Y29uamVjdHVyYWxfc2hhKEVsbGlwdGljQ3VydmUoSywgW2EgKyAxLCAwLCAxLCBhIC0gMiwgMF0pLCAxKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpUqLTHE/___code___.py", line 3, in <module>
    exec compile(u'conjectural_sha(EllipticCurve(K, [a + _sage_const_1 , _sage_const_0 , _sage_const_1 , a - _sage_const_2 , _sage_const_0 ]), _sage_const_1 )
  File "", line 1, in <module>
    
  File "/tmp/tmpx5SCfX/___code___.py", line 8, in conjectural_sha
    Q = saturate_prime(E, Q, ell)    
  File "/tmp/tmpoRYY4o/___code___.py", line 32, in saturate_prime
    Q2 = (ZZ(a2)*E(Q[_sage_const_0 ]) - ZZ(a1)*E(Q[_sage_const_1 ])).division_points(ell)[_sage_const_0 ]
IndexError: list index out of range