BSD + Zhang -- Sage

def EK_span(E, K, proof=False): ED = E.quadratic_twist(K.discriminant()) EK = E.change_ring(K) phi = ED.change_ring(K).isomorphism_to(EK) rat_gens = [EK(P) for P in E.gens(proof=proof)] twist_gens = [phi(P) for P in ED.gens(proof=proof)] return rat_gens + twist_gens def heegner_point_orbit(E, D, c, prec=200): print "Computing heegner point..." K.<a> = QuadraticField(D) y_c = E.heegner_point(D, c) P = y_c.point_exact(prec) print "Optimizing ring class field..." Kc, from_Kc, to_Kc = P[1].parent().optimized_representation() EKc = E.change_ring(Kc) P = EKc(to_Kc(P[0]), to_Kc(P[1])) if Kc.degree() < K.degree() * (c+1): print "Kc does not contain K, enlarging..." # Kc must contain K # There should be a better way to do this... # And the returned results are awful... Kc_rel.<b> = K.extension(Kc.polynomial()) Kc_new.<c> = Kc_rel.absolute_field() _, h_abs = Kc_new.structure() Kc_opt, _, h_opt = Kc_new.optimized_representation('a') Kc = Kc_opt EKc = E.change_ring(Kc) P = EKc((h_opt(h_abs(P[0])), h_opt(h_abs(P[1])))) print "K_%s = %s" % (c, Kc) print "y_%s = %s" % (c, P) K_to_Kc = K.hom([Kc(D).sqrt()]) EKc = E.change_ring(Kc) y_c = P a = Kc.gen() conjugates = Kc.polynomial().roots(Kc, multiplicities=False) # Dihedral extension over Q, want generator for cyclic part. for sigma_a in conjugates: sigma = Kc.hom([sigma_a]) if sigma(sigma(a)) == a: continue for p, e in factor(c+1): if (sigma**((c+1)//p))(a) == a: break else: good = sigma break # this is a generator return [EKc((sigma**k)(P[0]), (sigma**k)(P[1])) for k in range((c+1))] @cached_function def pair(P, Q): return ((P+Q).height() - P.height() - Q.height()) / 2

D = -7 K.<a> = QuadraticField(D)

E = EllipticCurve('389a') ED = E.quadratic_twist(D) EK = E.change_ring(K)

r = 3; c = 5; n = c+1

orbit = heegner_point_orbit(EllipticCurve('389a'), -7, 5)

Computing heegner point...
Optimizing ring class field...
K_5 = Number Field in a3 with defining polynomial x^12 - 4*x^11 + 6*x^10
- 5*x^9 + 5*x^8 - 9*x^7 + 21*x^6 - 9*x^5 + 5*x^4 - 5*x^3 + 6*x^2 - 4*x +
1
y_5 = (-19973/8505*a3^11 + 15728/1701*a3^10 - 22373/1701*a3^9 +
83717/8505*a3^8 - 29069/2835*a3^7 + 56684/2835*a3^6 - 130936/2835*a3^5 +
41431/2835*a3^4 - 39703/8505*a3^3 + 23101/1701*a3^2 - 17755/1701*a3 +
39217/8505 : 84071/59535*a3^11 - 53506/8505*a3^10 + 18131/1701*a3^9 -
563147/59535*a3^8 + 170729/19845*a3^7 - 313454/19845*a3^6 +
712174/19845*a3^5 - 485662/19845*a3^4 + 386977/59535*a3^3 -
617048/59535*a3^2 + 896264/59535*a3 - 534217/59535 : 1)

Computing heegner point...
Optimizing ring class field...
K_5 = Number Field in a3 with defining polynomial x^12 - 4*x^11 + 6*x^10 - 5*x^9 + 5*x^8 - 9*x^7 + 21*x^6 - 9*x^5 + 5*x^4 - 5*x^3 + 6*x^2 - 4*x + 1
y_5 = (-19973/8505*a3^11 + 15728/1701*a3^10 - 22373/1701*a3^9 + 83717/8505*a3^8 - 29069/2835*a3^7 + 56684/2835*a3^6 - 130936/2835*a3^5 + 41431/2835*a3^4 - 39703/8505*a3^3 + 23101/1701*a3^2 - 17755/1701*a3 + 39217/8505 : 84071/59535*a3^11 - 53506/8505*a3^10 + 18131/1701*a3^9 - 563147/59535*a3^8 + 170729/19845*a3^7 - 313454/19845*a3^6 + 712174/19845*a3^5 - 485662/19845*a3^4 + 386977/59535*a3^3 - 617048/59535*a3^2 + 896264/59535*a3 - 534217/59535 : 1)

EK_gens = EK_span(E, K)

h_mat = matrix(RR, 8, 8)

for i in range(3): for j in range(3): h_mat[i,j] = pair(EK_gens[i], EK_gens[j])

for i in range(c): for j in range(c): h_mat[i+3, j+3] = pair(orbit[0], orbit[j-i])

h_mat

[    0.686667083305587     0.268478098806726  6.66133814775094e-16    
0.000000000000000     0.000000000000000     0.000000000000000    
0.000000000000000     0.000000000000000]
[    0.268478098806726     0.327000773651605     0.000000000000000    
0.000000000000000     0.000000000000000     0.000000000000000    
0.000000000000000     0.000000000000000]
[ 6.10622663543836e-16 -2.77555756156289e-17      2.96766919125601    
0.000000000000000     0.000000000000000     0.000000000000000    
0.000000000000000     0.000000000000000]
[    0.000000000000000     0.000000000000000     0.000000000000000     
2.17273141179151     -1.07937746994395   -0.0538585294374920   
0.0937405869784393   -0.0538585294428808]
[    0.000000000000000     0.000000000000000     0.000000000000000    
-1.07937746994340      2.17273141179151     -1.07937746994395  
-0.0538585294374920    0.0937405869784393]
[    0.000000000000000     0.000000000000000     0.000000000000000  
-0.0538585294428808     -1.07937746994340      2.17273141179151    
-1.07937746994395   -0.0538585294374920]
[    0.000000000000000     0.000000000000000     0.000000000000000   
0.0937405869784393   -0.0538585294428808     -1.07937746994340     
2.17273141179151     -1.07937746994395]
[    0.000000000000000     0.000000000000000     0.000000000000000  
-0.0538585294374920    0.0937405869784393   -0.0538585294428808    
-1.07937746994340      2.17273141179151]

[    0.686667083305587     0.268478098806726  6.66133814775094e-16     0.000000000000000     0.000000000000000     0.000000000000000     0.000000000000000     0.000000000000000]
[    0.268478098806726     0.327000773651605     0.000000000000000     0.000000000000000     0.000000000000000     0.000000000000000     0.000000000000000     0.000000000000000]
[ 6.10622663543836e-16 -2.77555756156289e-17      2.96766919125601     0.000000000000000     0.000000000000000     0.000000000000000     0.000000000000000     0.000000000000000]
[    0.000000000000000     0.000000000000000     0.000000000000000      2.17273141179151     -1.07937746994395   -0.0538585294374920    0.0937405869784393   -0.0538585294428808]
[    0.000000000000000     0.000000000000000     0.000000000000000     -1.07937746994340      2.17273141179151     -1.07937746994395   -0.0538585294374920    0.0937405869784393]
[    0.000000000000000     0.000000000000000     0.000000000000000   -0.0538585294428808     -1.07937746994340      2.17273141179151     -1.07937746994395   -0.0538585294374920]
[    0.000000000000000     0.000000000000000     0.000000000000000    0.0937405869784393   -0.0538585294428808     -1.07937746994340      2.17273141179151     -1.07937746994395]
[    0.000000000000000     0.000000000000000     0.000000000000000   -0.0538585294374920    0.0937405869784393   -0.0538585294428808     -1.07937746994340      2.17273141179151]

Reg_Kc = h_mat.det(); Reg_Kc

3.99487194842550

3.99487194842550

Omega_Kc = (2*E.period_lattice().complex_area())^n; Omega_Kc

896795.322438472

896795.322438472

Delta_Kc = orbit[0][0].parent().disc(); Delta_Kc.factor()

5^10 * 7^6

5^10 * 7^6

lhs = factorial(r+c) * Omega_Kc * Reg_Kc / abs(Delta_Kc).sqrt().n(); lhs

134763.592281016

134763.592281016

leading_K = (E.lseries().taylor_series(1) * ED.lseries().taylor_series(1))[3]; leading_K

3.35867855873122

3.35867855873122

all = 2.9715e6 # product of L'(f, chi, 1), computed below

rhs = leading_K * all; rhs

9.98031333726981e6

9.98031333726981e6

rhs / lhs

74.0579348497799

74.0579348497799

E = EllipticCurve('389a') yc = E.heegner_point(-7,5) Kc = yc.ring_class_field() K = Kc.quadratic_field() G = Kc.galois_group(K); G # allQ = [g.quadratic_form() for g in G] if G.cardinality() > 1: allQ = [prod([G[1]]*k, G[0]).quadratic_form() for k in range(G.cardinality())] else: allQ = G[0].quadratic_form() Q = allQ[0] N = E.conductor(); D = Q.discriminant()

Lgz = Dokchitser(N**2*D**2, [0,0,1,1], weight=2, eps=-1, prec=20) nterms = ZZ(Lgz.gp()('cflength()')); print nterms

from sage.modular.modform.l_series_coeffs import gross_zagier_L_series nterms = 420000 an_list = E.anlist(nterms) L_A_list = [gross_zagier_L_series(an_list, Q, N, 2, 'q') for Q in allQ] C.<zeta> = CyclotomicField(c+1)

L_A_list = [L_A.change_ring(C) for L_A in L_A_list]

def L_chi_prod(): all = 1 for i in [1..5]: print "i", i coeffs = list(sum(zeta^(i*k)*L_A for k, L_A in enumerate(L_A_list)).change_ring(ZZ)) print "got coeffs..." Lgz.gp().set('a', coeffs[1:]) Lgz.init_coeffs('a[k]', 1) print "check", Lgz.check_functional_equation() val = Lgz.derivative(1) print "val", val all *= val return all

all = L_chi_prod()

i 1
got coeffs...
check 5.8673e-6
val 9.3843
i 2
got coeffs...
check 7.1526e-7
val 30.284
i 3
got coeffs...
check 0.000010490
val 36.790
i 4
got coeffs...
check 7.1526e-7
val 30.284
i 5
got coeffs...
check 5.8673e-6
val 9.3843

i 1
got coeffs...
check 5.8673e-6
val 9.3843
i 2
got coeffs...
check 7.1526e-7
val 30.284
i 3
got coeffs...
check 0.000010490
val 36.790
i 4
got coeffs...
check 7.1526e-7
val 30.284
i 5
got coeffs...
check 5.8673e-6
val 9.3843

all

2.9715e6

2.9715e6

all = 2.9715e6

rhs = leading_K * all; rhs

9.98031333726981e6

9.98031333726981e6

rhs / lhs

4739.70783038591

4739.70783038591

KcP = orbit[0][0].parent() for p in primes(25): print p print E.division_polynomial(p).change_ring(QQ).roots(KcP)

2
[]
3
[]
5
[]
7
[]
11
[]
13
[]
17
[]
19
[]
23
[]

2
[]
3
[]
5
[]
7
[]
11
[]
13
[]
17
[]
19
[]
23
[]

rhs = leading_K * L_chi_prod()

i 1
got coeffs...
check 5.566e-6
val 9.384
i 2
got coeffs...
check 0.00003195
val 30.28
i 3
got coeffs...
check -0.00004005
val 36.79
i 4
got coeffs...
check 0.00003195
val 30.28
i 5
got coeffs...
check 5.566e-6
val 9.384

i 1
got coeffs...
check 5.566e-6
val 9.384
i 2
got coeffs...
check 0.00003195
val 30.28
i 3
got coeffs...
check -0.00004005
val 36.79
i 4
got coeffs...
check 0.00003195
val 30.28
i 5
got coeffs...
check 5.566e-6
val 9.384

rhs = leading_K * 9.384^2 * 30.28^2 * 36.79; rhs

9.97667832263972e6

9.97667832263972e6

9.97667832263972e6 / 98734.3766370173

101.045640459326

101.045640459326

%time zeta * L_A_list[0].change_ring(C)

^C
^C
Traceback (click to the left of this block for traceback)
...
__SAGE__

^C
^C
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_138.py", line 11, in <module>
    open("___code___.py","w").write("# -*- coding: utf-8 -*-\n" + _support_.preparse_worksheet_cell(base64.b64decode("emV0YSAqIExfQV9saXN0WzBd"),globals())+"\n"); execfile(os.path.abspath("___code___.py")); print "CPU time: %.2f s,  Wall time: %.2f s"%(cputime(__SAGE_t__), walltime(__SAGE_w__))
  File "", line 1, in <module>
    
  File "/private/var/folders/o3/o3m2i3FKHHS3wWaoXqhwYk+++TU/-Tmp-/tmp9ljUxE/___code___.py", line 3, in <module>
    zeta * L_A_list[_sage_const_0 ]
  File "", line 1, in <module>
    
  File "/Users/robertwb/sage/sage-4.0/local/lib/python2.6/site-packages/sage/misc/displayhook.py", line 174, in displayhook
    print_obj(sys.stdout, obj)
  File "/Users/robertwb/sage/sage-4.0/local/lib/python2.6/site-packages/sage/misc/displayhook.py", line 142, in print_obj
    print >>out_stream, `obj`
  File "sage_object.pyx", line 101, in sage.structure.sage_object.SageObject.__repr__ (sage/structure/sage_object.c:1370)
  File "power_series_ring_element.pyx", line 502, in sage.rings.power_series_ring_element.PowerSeries._repr_ (sage/rings/power_series_ring_element.c:4188)
  File "number_field_element.pyx", line 438, in sage.rings.number_field.number_field_element.NumberFieldElement.__repr__ (sage/rings/number_field/number_field_element.cpp:6257)
  File "number_field_element.pyx", line 1738, in sage.rings.number_field.number_field_element.NumberFieldElement.polynomial (sage/rings/number_field/number_field_element.cpp:13331)
  File "parent.pyx", line 811, in sage.structure.parent.Parent.__call__ (sage/structure/parent.c:6188)
  File "coerce_maps.pyx", line 82, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (sage/structure/coerce_maps.c:3136)
  File "coerce_maps.pyx", line 77, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (sage/structure/coerce_maps.c:3038)
  File "/Users/robertwb/sage/sage-4.0/local/lib/python2.6/site-packages/sage/rings/polynomial/polynomial_ring.py", line 312, in _element_constructor_
    return C(self, x, check, is_gen, construct=construct, **kwds)
  File "/Users/robertwb/sage/sage-4.0/local/lib/python2.6/site-packages/sage/rings/polynomial/polynomial_element_generic.py", line 664, in __init__
    self.__poly = pari(str(x)).Polrev()
KeyboardInterrupt
__SAGE__

%time v = L_A_list[0].change_ring(C)

CPU time: 12.58 s,  Wall time: 12.80 s

CPU time: 12.58 s,  Wall time: 12.80 s

%time v = zeta * v

CPU time: 1.58 s,  Wall time: 1.62 s

CPU time: 1.58 s,  Wall time: 1.62 s

chi1 = list(sum(zeta^(k)*L_A for k, L_A in enumerate(L_A_list)).change_ring(ZZ))

import cPickle as pickle pickle.dump(chi1, open("chi_1.sobj", 'w'))

chi_1.sobj

def gz_via_lcalc(an_list): f = open("data.txt", "w") f.write("2\n") # real coeffs f.write("0\n") # unspecified L-func f.write("%s\n" % (len(an_list)-1)) # number of coeffs f.write("0\n") # non-periodic coeffs f.write("2\n") # gamma factor count f.write("1\n0 0\n1\n0 0\n") # Gamma(s)^2 f.write("%s\n" % RR(D*N/(4*pi)^2)) # Q f.write("-1 0\n") # epsilon f.write("0\n") # poles for n,a in enumerate(an_list): # must center critical line at 1/2 if n == 0: continue f.write("%s\n" % (a*sqrt(RR(n)))) f.close() res = lcalc("-v -x .5 -y 0 -F data.txt") print res return RR(res)

gz_via_lcalc([1] * 100)

lcalc(17638) malloc: *** mmap(size=4292071424) failed (error code=12)
*** error: can't allocate region
*** set a breakpoint in malloc_error_break to debug
terminate called after throwing an instance of 'std::bad_alloc'
  what():  St9bad_alloc

Traceback (click to the left of this block for traceback)
...
TypeError: Unable to convert x (='') to real number.

lcalc(17638) malloc: *** mmap(size=4292071424) failed (error code=12)
*** error: can't allocate region
*** set a breakpoint in malloc_error_break to debug
terminate called after throwing an instance of 'std::bad_alloc'
  what():  St9bad_alloc

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_202.py", line 9, in <module>
    open("___code___.py","w").write("# -*- coding: utf-8 -*-\n" + _support_.preparse_worksheet_cell(base64.b64decode("Z3pfdmlhX2xjYWxjKFsxXSAqIDEwMCk="),globals())+"\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/private/var/folders/o3/o3m2i3FKHHS3wWaoXqhwYk+++TU/-Tmp-/tmpQkP4o1/___code___.py", line 3, in <module>
    gz_via_lcalc([_sage_const_1 ] * _sage_const_100 )
  File "", line 1, in <module>
    
  File "/private/var/folders/o3/o3m2i3FKHHS3wWaoXqhwYk+++TU/-Tmp-/tmp9252LF/___code___.py", line 22, in gz_via_lcalc
    return RR(res)
  File "parent.pyx", line 811, in sage.structure.parent.Parent.__call__ (sage/structure/parent.c:6188)
  File "coerce_maps.pyx", line 82, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (sage/structure/coerce_maps.c:3136)
  File "coerce_maps.pyx", line 77, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (sage/structure/coerce_maps.c:3038)
  File "real_mpfr.pyx", line 404, in sage.rings.real_mpfr.RealField_class._element_constructor_ (sage/rings/real_mpfr.c:5075)
  File "real_mpfr.pyx", line 1027, in sage.rings.real_mpfr.RealNumber._set (sage/rings/real_mpfr.c:8816)
TypeError: Unable to convert x (='') to real number.

data.txt

BSD + Zhang

655 days ago by robertwb