rank_two_computation_sha_reg_MW

K.<a> = NumberField(x^2 - x - 1) import psage.ellcurve.lseries.sqrt5 as sqrt5w 

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

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(40): 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 

time conjectural_sha(E, 1.44144111075908) 

2
3
5
7
11
13
17
19
23
29
31
37
1.00000000000001
Time: CPU 69.62 s, Wall: 69.63 s

2
3
5
7
11
13
17
19
23
29
31
37
1.00000000000001
Time: CPU 69.62 s, Wall: 69.63 s