Dembele's algorithm (for modular elliptic curves over real quadratic fields)

var('x') K.<a> = NumberField(x^2 - x - 1) R = RealField(200) C = ComplexField(200) phi2, phi = K.embeddings(R) # I think it's important that we choose the correct embedding here. # We want the embedding such that phi(a) > 0 class QuadraticCharacter: def __init__(self, P): P = ideal(P) self.k = P.residue_field() self.K = self.k.lift(self.k(1)).parent().number_field() self.character_table = {} self.character_table[self.k(0)] = 0 self.conductor = P g = self.k.multiplicative_generator() a = g current_sign = -1 for n in range(self.k.order() - 1): self.character_table[a] = current_sign current_sign = -current_sign a = a * g # a should now be gen. We use this to make sure we haven't made a big mistake: if a != g: print a, g raise ValueError("something went wrong.") def __call__(self, a): a = self.k(a) return self.character_table[a] def gauss_sum(self): CC = ComplexField(200) RR = RealField(200) I = CC.0 characteristic = self.k.characteristic() x = CC(0) d = K(K.disc()).sqrt() m = self.conductor.gens_reduced()[0] for a, value in self.character_table.iteritems(): x = x + value * CC(2 * pi * I * (self.k.lift(a)/(d * m)).trace()).exp() return x def make_tot_pos(n): if phi(n) > 0: if phi2(n) < 0: return n * a else: return n else: if phi2(n) < 0: return n * a * (1 - a) else: return n * (1 - a) 

E = EllipticCurve(K, [1, -(1-a + 1), 1-a, 0, 0]) # this is the curve we are going to be looking for # We are now about to cheat, and use the curve that we are looking for to compute a bunch of ap values ap_list = {} ideals = K.ideals_of_bdd_norm(2000) for n in range(1, 2001): for p in ideals[n]: p = K.ideal(p.gens_reduced()) if p == E.conductor(): ap_list[p] = -1 elif is_prime(p): ap_list[p] = E.change_ring(p.residue_field()).trace_of_frobenius() an_list = {} for n in range(1, 2001): for I in ideals[n]: an = 1 F = I.factor() for (p, e) in F: p = K.ideal(p.gens_reduced()) ap = ap_list[p] if p == E.conductor(): an *= ap^e else: if e == 1: an *= ap elif e == 2: an *= (ap^2 - norm(p)) elif e == 3: an *= (ap^3 - 2 * norm(p) * ap) elif e == 4: an *= (ap^4 - 4 * ap^2 * norm(p) + norm(p)^2) elif e == 5: an *= ap^5 - 4 * ap^3 * norm(p) + 3 * ap * norm(p)^2 else: raise ValueError("not ready to handle thing that are divisible by higher than cubic powers") an_list[I] = an 

def mixed_period(f, chi, normbound = 10000000): # right now f is the form with the above ap list N = phi(a * (5 * a - 2)) # a totally positive generator for the the level of f Nbar = phi2(a * (5 * a - 2)) c = make_tot_pos(chi.conductor.gens_reduced()[0]) nu = phi(c) nubar = phi2(c) D = R(5) S = 0 for (n, an) in an_list.iteritems(): if norm(n) > normbound: continue n = make_tot_pos(n.gens_reduced()[0]) S2 = 0 for k in range(-10, 10): # not sure what kind of range to use on k right now A = R(-2 * pi * phi(n) * phi(a)^(2 * k + 1) / (nu * sqrt(D * N))) B = R(2 * pi / sqrt(D) * (phi2(n) * phi(1 - a)^(2 * k + 1) / (nubar * sqrt(Nbar)) - phi(n) * phi(a)^(2 * k)/(nu * sqrt(N)))) S2 = S2 + (1 - A.exp())*B.exp() S = S + an * chi(n)/norm(n) * S2 S = -S * 5/R(2 * pi^2) return C(-4 * pi^2 * sqrt(5) * chi.gauss_sum()) * S 

# now we define some characters: chi_mm = QuadraticCharacter(K(3)) chi_mp = QuadraticCharacter(K.primes_above(11)[0]) chi_pm = QuadraticCharacter(K.primes_above(11)[1]) # the names of the characters come from their values on a and 1 - a: print chi_mm(a), chi_mm(1 - a) print chi_mp(a), chi_mp(1 - a) print chi_pm(a), chi_pm(1 - a) 

-1 -1
-1 1
1 -1

-1 -1
-1 1
1 -1

# now we compute some mixed periods with these characters mm = mixed_period(None, chi_mm, 1000); print mm/5 mp = mixed_period(None, chi_mp, 1000); print mp/5 # notice that we happened to guess that we should divide by 2 pm = mixed_period(None, chi_pm, 1000)/2; print pm/5 # we divide by 5 when printing to match the values with Lassina's (don't know right now which values are correct, but it doesn't matter) 

7.5426808112369306740153992460705909078327594653694558750026 +
3.1292145282904993254146053918433466137978573813023609226541e-60*I
5.2221350999144463281690248456963896493064281527625362940368e-59 +
20.241415522709212230291917129318279518266661747418781293909*I
2.4759420110586928347012906323886284823668688899932755004223e-59 +
9.5969503080822906114254153871418851015705865742327501285497*I

7.5426808112369306740153992460705909078327594653694558750026 + 3.1292145282904993254146053918433466137978573813023609226541e-60*I
5.2221350999144463281690248456963896493064281527625362940368e-59 + 20.241415522709212230291917129318279518266661747418781293909*I
2.4759420110586928347012906323886284823668688899932755004223e-59 + 9.5969503080822906114254153871418851015705865742327501285497*I

# the other mixed period can be easily computed from those 3: pp = -mp * pm / mm; print pp/5 

25.754219726133373182316881287139936136719897307889005011941 -
1.4357254757988874927592999495936713823118020704451155727731e-58*I

25.754219726133373182316881287139936136719897307889005011941 - 1.4357254757988874927592999495936713823118020704451155727731e-58*I

tau1 = -pp/mp; print tau1 tau2 = 1/2 * (1 - pp/pm); print tau2 

3.8104336361382328973946083064271792914951703343080270986463e-60 +
1.2723527016793487415983952188444507211647829380853493895806*I
0.50000000000000000000000000000000000000000000000000000000000 +
1.3417918661329250385953663654972825159308341044426131079028*I

3.8104336361382328973946083064271792914951703343080270986463e-60 + 1.2723527016793487415983952188444507211647829380853493895806*I
0.50000000000000000000000000000000000000000000000000000000000 + 1.3417918661329250385953663654972825159308341044426131079028*I

# now compute the j-invariants j1 = j_invariant_qexp(10)( CC(2 * pi * i * tau1).exp() ).real(); print j1 j2 = j_invariant_qexp(10)( CC(2 * pi * i * tau2).exp() ).real(); print j2 # there are some small differences from Lassina's values, we are slightly closer to the actual values than he is, though. 

3777.28932046422
-3883.75245635613

3777.28932046422
-3883.75245635613

print phi2(E.j_invariant()) - j1, phi(E.j_invariant()) - j2 #something funny might be going on with which j invariant is which 

-0.0262921690950861 0.102331286809203

-0.0262921690950861 0.102331286809203

D = E.discriminant(); print D 

-16*a + 25

-16*a + 25

# now we compute c4a = (j1 * phi2(D))^(1/3) c4b = (j2 * phi(D))^(1/3) print c4a, c4b 

50.8886618920066 15.1115889040392

50.8886618920066 15.1115889040392

# now we happen to notice that if we take c4 = 41 - 16 * a #then print phi2(c4) - c4a, phi(c4) - c4b #these values are almost 0, so this must be the right value of c4 

-0.000118072008277181 -0.000132724037516851

-0.000118072008277181 -0.000132724037516851

#now compute c6 = sqrt(c4^3 - 1728 * D) 

# and finally, take E2 = EllipticCurve_from_c4c6(c4, c6) 

# to check, we see that E2 is the curve that we started from: print E.is_isomorphic(E2) 

True

True