#################################################################################
#
# (c) Copyright 2011 William Stein
#
# This file is part of PSAGE
#
# PSAGE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# PSAGE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see .
#
#################################################################################
"""
Computation of Chow-Heegner points using Zhang's construction. This
is for an appendix by William Stein to a paper by
Darmon-Daub-Lichtenstein-Rotger.
TODO:
[ ] when finding representatives in upper half plane, always take
one with largest imaginary part
"""
from sage.all import *
def chow_heegner1(E, F, num=1, prec=53, ser_prec=50, alpha=None):
K = ComplexField(prec) if prec > 53 else CDF
q = K['q'].gen()
h = sum(F.an(n)/n * q**n for n in range(1,ser_prec+1))
if alpha is not None:
h -= alpha
# Roots in the unit disk
roots = [z for z,_ in h.roots(K) if abs(z) < 1]
# Roots in the plane
pi = K.pi(); I = K.gen()
roots = [z.log()/(2*pi*I) for z in roots if z]
# The ones in the upper half plane
roots = [z for z in roots if z.imag() > 0]
# Sort by imaginary part (largest first)
roots.sort(lambda a,b: cmp(b[1],a[1]))
# Map down
phi_E = E.modular_parametrization()
return [(r, phi_E(r)) for r in roots[:num]]
def multiple_of_gen(E, P):
prec = 53
Q = E.gens()[0]
while True:
s = sqrt(P.height(prec) / Q.height(prec))
try:
m = ZZ(QQ(s))
if (m * Q - P).order() > 0:
return m
except TypeError: pass
prec += 10
if prec >= 1000:
raise RuntimeError
def recognize(E, P, n, prec=53, **kwds):
v = []
rel = algebraic_dependency(P[0], n, **kwds)
print "rel = ", rel
if n % rel.degree() != 0:
print "WARNING: n = %s, but degree = %s"%(n, rel.degree())
for f, k in rel.factor():
g, d = sage.schemes.elliptic_curves.heegner.make_monic(f)
K = NumberField(g, names='a')
a = K.gen()
E1 = E.change_ring(K)
lifts = E1.lift_x(a/d, all=True)
if len(lifts) == 0:
raise ValueError, "no point over field defined by %s with X-coordinate %s"%(
g, a/d)
C = ComplexField(prec) if prec > 53 else CDF
phi = K.embeddings(C)[0]
v = [(abs(phi(lifts[i][1]) - P[1]), i) for i in range(len(lifts))]
v.sort()
Q = lifts[v[0][1]]
E2 = E.change_ring(CDF)
R = sum( [E2.point([e(Q[0]),e(Q[1]),e(Q[2])], check=False) for e in K.embeddings(C)] )
prec2 = prec
R0 = 'no point found'
while prec2>0:
M = RealField(prec2)
try:
R0 = E([M(z[0]) for z in list(R)])
break
except:
prec2 -= 1
m = multiple_of_gen(E, R0)
v.append((k, m, n/rel.degree(), R0, "(%s)^%s"%(f,k), R, f, Q))
return v
######################################################
from sage.schemes.elliptic_curves.heegner import make_monic
from sage.rings.all import (QQ, ZZ, CDF, RDF,
NumberField, ComplexField, is_ComplexField)
def lift_x_rational(E, x, var='a'):
"""
Given a rational number x, return list of the points on E, over
either the rational numbers or a quadratic extension, with that x
coordinate.
EXAMPLES::
sage: E = EllipticCurve([1..5])
sage: from psage.ellcurve.points.chow_heegner import lift_x
sage: lift_x_rational(E, 1)
[(1 : 2 : 1), (1 : -6 : 1)]
sage: v = lift_x_rational(E, 2); v
[(2 : a : 1), (2 : -a - 5 : 1)]
sage: sum(v)
(0 : 1 : 0)
"""
R = QQ['y']
f = E.defining_polynomial()
g = f.parent().hom([x, R.gen(), 1])(f)
if g.is_irreducible():
# point is over quadratic extension
h, d = make_monic(g)
K = NumberField(h, names=var)
a = K.gen()
y = a/d # one solution
P = E.change_ring(K)([x,y])
Q = -P
v = [P]
if Q != P:
v.append(Q)
return v
else:
# point is on curve
return [E([x,y]) for y in g.roots(multiplicities=False)]
def lift_y(E, y, var='a'):
"""
Given a number y, return list all of the points on E over the base
field of E with that y-coordinate.
EXAMPLES::
"""
K = E.base_ring()
R = K['x']
f = E.defining_polynomial()
g = f.parent().hom([R.gen(), y, 1])(f)
if g.is_irreducible():
# point is over quadratic extension
return []
else:
# point is on curve over base
return [E([x,y]) for x in g.roots(multiplicities=False)]
def elliptic_logarithms(P, prec=53):
"""
Return the elliptic logarithms of P to prec bits of precision, for
each of the complex embeddings of the base field.
INPUT:
- P -- point on elliptic curve over QQ or number field
- prec -- integer; bits of precision
"""
E = P.curve()
return [E.period_lattice(phi).elliptic_logarithm(P, prec=prec)
for phi in E.base_ring().embeddings(ComplexField(prec))]
def B_bound(ymin, prec):
"""
Return integer B so that using B terms of the series of a modular
parameterization the tail end of the sum is 0 for any point with y
at least ymin.
INPUT:
- ``ymin`` -- positive real number; minimum y coordinate
- ``prec`` -- positive integer (bits of precision)
"""
y = RR(ymin)
epsilon = RR(2)**(-(prec+1))
pi = RR.pi()
return int((epsilon*(1 - (-2*pi*y).exp())).log() / (-2*pi*y)) + 1
def refine_root(f, root, max_steps=100, verbose=False):
"""
Use the given number of steps of the Newton-Raphson algorithm
to refine the given input root of f.
INPUT:
- `f` -- a polynomial
- ``root`` -- a complex number
- ``steps`` -- a positive integer
"""
C = f.base_ring()
tiny = C(2)**(5 - C.prec())
f_prime = f.derivative()
#if f.degree() < 25000: # see trac 11766: bug in fast_callable. :-(
# f = fast_callable(f, root.parent())
# f_prime = fast_callable(f_prime, root.parent())
if 1:
f = ComplexPolynomial(f)
f_prime = ComplexPolynomial(f_prime)
last_root = root
for i in range(max_steps):
root = root - f(root)/f_prime(root)
if abs(last_root - root) < tiny:
print "early break after %s iter: f(root) = %s"%(i, f(root))
break
else:
last_root = root
return root
def exact_lifts(E, f):
"""
INPUT:
- `E` -- elliptic curve over QQ
- `f` -- polynomial that x(Q) approximately satisfies
OUTPUT:
- the 0, 1 or 2 points on E over number field with x
coordinate a root of f.
"""
g, d = sage.schemes.elliptic_curves.heegner.make_monic(f)
K = NumberField(g, assume_disc_small=True, names='a')
a = K.gen()
return E.change_ring(K).lift_x(a/d, all=True)
def exact_point(E, Q, f):
"""
INPUT:
- `E` -- elliptic curve over QQ
- `Q` -- point on base extension of E over some complex
floating point field
- `f` -- polynomial that x(Q) approximately satisfies
OUTPUT:
- point on E over number field that maps to Q under
one of the embeddings, or ValueError, if none exists
"""
v = exact_lifts(E, f)
if len(v) == 0:
raise ValueError, "there is no point with given x coordinate over field, so f=%s is wrong"%f
if len(v) == 1:
return v[0]
K = v[0].curve().base_field()
# Find psi such that psi(v[0]) = x(Q)
e = K.embeddings(Q.curve().base_field())
if len(e) == 1 or (abs(e[0](v[0][0]) - Q[0]) < abs(e[1](v[0][0]) - Q[0])):
psi = e[0]
else:
psi = e[1]
if abs(psi(v[0][1]) - Q[1]) < abs(psi(v[1][1]) - Q[1]):
R = v[0]
else:
R = v[1]
return R
def exact_lifts_y(E, f):
g, d = sage.schemes.elliptic_curves.heegner.make_monic(f)
K = NumberField(g, assume_disc_small=True, names='a')
a = K.gen()
return lift_y(E.change_ring(K), a/d)
def exact_point_y(E, Q, f):
v = exact_lifts_y(E, f)
if len(v) == 0:
raise ValueError, "there is no point with given y coordinate over field, so f=%s is wrong"%f
if len(v) == 1:
return v[0]
K = v[0].curve().base_field()
# Find psi such that psi(v[0]) = x(Q)
e = K.embeddings(Q.curve().base_field())
if len(e) == 1 or (abs(e[0](v[0][1]) - Q[1]) < abs(e[1](v[0][1]) - Q[1])):
psi = e[0]
else:
psi = e[1]
if abs(psi(v[0][0]) - Q[0]) < abs(psi(v[1][0]) - Q[0]):
R = v[0]
else:
R = v[1]
return R
from psage.ellcurve.points.chow_heegner_fast import (
cdf_roots_of_rdf_poly, Polynomial_RDF_gsl, ComplexPolynomial)
def zeros_in_h(f):
K = f.base_ring()
# Roots in the unit disk
if K.prec() > 53:
zeros = [z for z,_ in f.roots(K) if abs(z) < 1]
else:
# find roots to double precision, then coerce into K
zeros = [z for z in f.roots(CDF,multiplicities=False) if abs(z) < 1]
if K != CDF:
zeros = [K(z) for z in zeros]
# Zeros in the plane
pi = K.pi(); I = K.gen()
zeros = [z.log()/(2*pi*I) for z in zeros if z]
# The ones in the upper half plane
zeros = [z for z in zeros if z.imag() > 0]
# Sort by imaginary part (largest first)
zeros.sort(lambda a,b: cmp(b[1],a[1]))
return zeros
def point_in_fiber(P, prec, index=0, embedding=0, low_series_prec=50, verbose=False):
"""
INPUT:
- P -- point on elliptic curve E over a number field
- prec -- bits of precision
- index -- which low precision point to start with in upper
half plane; index=0 is the highest point, index=1 the next
highest, etc., ordered by y coordinate.
- embedding -- which embedding of number field
into complex field to use
- low_series_prec -- optional low series precision used to
find low precision point in fiber.
OUTPUT:
- point z in the upper half plane to precision prec that
maps to P.
"""
E = P.curve().change_ring(QQ)
# TODO: Often elliptic logarithm totally hangs; using higher precision fixes that. :-(
# SEE trac 11767: http://trac.sagemath.org/sage_trac/ticket/11767
if verbose: print "elliptic logarithms..."
#print "starting elliptic logarithm"
#print E, P, 100+prec
P0 = elliptic_logarithms(P, 100+prec)[embedding]
#print "didn't hang"
# Find z to very low precision
f0 = phi_poly(E, low_series_prec, base_field=CDF)
if verbose: print "low precision degree = ", f0.degree()
z = zeros_in_h(f0 - P0)
for i, x in enumerate(z):
print i, x.imag()
try:
tau0 = z[index]
except IndexError:
raise ValueError, "no point with index %s to given precision in upper half plane (not enough points)"%index
if tau0.imag() <= 0:
raise ValueError, "no point with index %s not in upper half plane"%index
# Refine z to be a root to the requested precision
C = ComplexField(prec)
f1 = phi_poly(E, B_bound(tau0.imag(), prec=prec), base_field=C)
pi = CDF.pi(); I = CDF.gen()
q0 = CDF(exp(2*pi*I*tau0))
if True or verbose: print "refinement degree = %s (imag = %s)"%(f1.degree(), tau0.imag())
q1 = refine_root(f1 - P0, C(q0), verbose=verbose)
if verbose: print "refined."
tau1 = q1.log()/(2*C.pi()*C.gen())
return tau1
def phi_poly(E, B, base_field=QQ):
R = base_field['q']
v = E.anlist(B+1)
return R([0] + [v[n]/n for n in range(1,B+1)])
def absolute_trace_numerical(P, prec):
"""
Numerically compute the trace of the point P over an absolute
number field, using the given numbers of bits of precision.
INPUT:
- P -- point over a number field on an elliptic curve
that is actually defined over QQ.
- prec -- positive integer; bits of precision
"""
C = ComplexField(prec)
E = P.curve()
K = E.base_ring()
E2 = EllipticCurve([C(QQ(a)) for a in E.a_invariants()])
R = sum( [E2.point([e(P[0]),e(P[1]),e(P[2])], check=False)
for e in K.embeddings(C)] )
return R
def frob(P, i):
if i == 0: return P
E = P.curve()
k = E.base_field()
p = k.characteristic()
q = p**i
return E([a**q for a in list(P)])
def next_good_prime(D, p):
p = next_prime(p)
while D%p == 0:
p = next_prime(p)
return p
def absolute_trace_multimodular(P, stabilize=6):
E = P.curve()
EQQ = E.change_ring(QQ)
K = E.base_field()
D = ZZ(K.defining_polynomial().discriminant())
p = next_good_prime(D, 10007)
C = [Mod(0,1), Mod(0,1)] # so far
L = None # last
i = 0
while True:
S = 0
Ep = E.change_ring(GF(p))
for r in K.primes_above(p):
k = r.residue_field()
Pmod = E.change_ring(k)(P)
S += Ep(sum(frob(Pmod,i) for i in range(k.degree())))
# Now use CRT to refine our guess C:
if S[2] != 0: # 0 point
R = IntegerModRing(p)
C[0] = C[0].crt(R(S[0]))
C[1] = C[1].crt(R(S[1]))
p = next_good_prime(D, p)
i += 1
# test to see if it stabilized
if i % stabilize == 0:
if C[0].modulus() == 1 and C[1].modulus() == 1:
# always got 0, so point is 0.
return EQQ(0)
try:
L_new = [C[0].rational_reconstruction(), C[1].rational_reconstruction()]
if L is not None:
if L == L_new:
return EQQ(L_new)
L = L_new
except ValueError, msg:
pass
def label(E):
try:
return E.cremona_label()
except RuntimeError:
return '%s?'%(E.conductor(),)
class ChowHeegnerPoint(object):
def __init__(self, E, F):
if E.conductor() != F.conductor():
raise NotImplementedError
# curves should be optimal -- some checking.
assert label(E)[-1] in ['?','1']
assert label(F)[-1] in ['?','1']
self.E = E
self.F = F
def __repr__(self):
return "Chow-Heegner point on %s associated to %s"%(
label(self.E), label(self.F))
def points(self, i=4, starting_prec=500, starting_low_series_prec=150,
max_prec=1000, verbose=False, verbose1=True):
E = self.E; F = self.F
if verbose1:
print "E=%s, F=%s, deg(phi_F)=%s"%(
E.cremona_label(), F.cremona_label(), F.modular_degree())
v = []
m = self.F.modular_degree()
xP = 0
prec = starting_prec
low_series_prec = starting_low_series_prec
while i >= 0:
try:
P, f, z = self.point0(xP=xP, prec=prec, low_series_prec=low_series_prec, verbose=verbose)
k = m / f.degree()
P = k.numerator() * (P.division_points(k.denominator())[0])
n = round((P.height() / self.E.regulator()).sqrt())
data = {'xP':xP, 'index':n,
'P':P, 'f':f,
'prec':prec, 'low_series_prec':low_series_prec,
'degF':m, 'status':'success'}
if verbose1 or verbose: print data
v.append(data)
update = True
except ValueError, msg:
print msg
if prec+500 > max_prec:
data = {'xP':xP,'prec':prec, 'low_series_prec':low_series_prec,
'degF':m, 'status':'fail'}
if verbose1 or verbose: print data
v.append(data)
update = True
else:
prec += 500
low_series_prec += 30
print "Retry with prec=%s, low_series_prec=%s"%(prec, low_series_prec)
update = False
if update:
i -= 1
if xP > 0:
xP = -xP
else:
xP = -xP + 1
prec = starting_prec
low_series_prec = starting_low_series_prec
return v
def point0(self, xP=0, index=0, prec=500, low_series_prec=200, verbose=False):
"""
Compute (not provably) the Chow-Heegner point as a rational
point on E using the following algorithm.
1. Compute point P=(x,y) in F, starting with x=1, but
increasing as necessary.
2. Find point z to prec bits on X_0(N) that maps to P.
3. Map z to point Q in E.
4. Use LLL to find the best polynomial f of degree
m = deg(phi_F) * [QQ(P):QQ]
that is satisfied by x(Q). If f is not
irreducible, choose a different random point P in 1.
5. Compute the two points on E with x-coordinate x(Q),
and let R be the one that is numerically close to Q.
If there are no such points, increase prec and try again.
6. Compute the trace of R down to the rational numbers using
a numerical or multimodular algorithm. Output this trace.
"""
v = verbose
E = self.E; F = self.F
d = F.modular_degree()
while True:
# step 1
if v: print "step 1: compute P on F"; t=walltime()
P = lift_x_rational(F, xP)[0]
if v: print "P =", P
m = d * P.curve().base_field().degree()
if v: print walltime(t)
# step 2
if v: print "step 2: Find point z to %s bits prec that maps to P"%prec; t=walltime()
z = point_in_fiber(P, prec, index, low_series_prec=low_series_prec, verbose=v)
if v: print "z = ", z
if v: print walltime(t)
# step 3
if v: print "step 3: Map z to a point Q in E using the Weierstrass P function"; t=walltime()
# This step may dominate timewise.
Q = E.modular_parametrization()(z)
if v: print walltime(t)
try:
# step 4
if v: print "step 4: Find best polynomial satisfied by x(Q)"; t=walltime()
f = algebraic_dependency(Q[0], m, known_bits=prec)
# Take largest degree factor
f = f.factor()[-1][0]
if v: print "f =", f
if v: print "f(Q[0]) =", f(Q[0])
if v: print walltime(t)
# step 5
if v: print "step 5: compute best exact point R on E corresponding to Q"; t=walltime()
R = exact_point(E, Q, f)
if v: print walltime(t)
except ValueError:
# Try y-coordinate instead
# step 4
if v: print "step 4b: Find best polynomial satisfied by y(Q)"; t=walltime()
f = algebraic_dependency(Q[1], m, known_bits=prec)
# Take largest degree factor
f = f.factor()[-1][0]
if v: print "f =", f
if v: print "f(Q[1]) =", f(Q[1])
if v: print walltime(t)
# step 5
if v: print "step 5b: compute best exact point R on E corresponding to Q"; t=walltime()
R = exact_point_y(E, Q, f)
if v: print walltime(t)
# step 6
if v: print "step 6: trace R to E(QQ)"; t=walltime()
T = absolute_trace_multimodular(R)
if v: print walltime(t)
return T, f, z
#N = absolute_trace_numerical(R, prec)
# divide T by 2
if f.degree() != d:
T2 = T.division_points(f.degree() // d)
N2 = N.division_points(f.degree() // d)
if len(T2) == 1:
return T2[0]
else:
return T2, ('WARNING', N2, R)
# TODO: In this case, resolve the ambiguity numerically.
raise NotImplementedError
#return point_trace_numerical(R, prec-10, recognize=False)
"""
Next thing to try!
Since we are randomizing our input point P, we can assume that the whole
fiber phi_F^(-1)(P) maps to deg(phi_F) *distinct* points on E.
Thus if we compute a few points in the fiber to sufficient precision,
we will know for a fact we've found the whole fiber. This way we
can hopefully deal with 153a1, 153b1.
"""
######################################################################################
# Gamma0(N) equivalence of points in the upper half plane
######################################################################################
def sl2z_rep_in_fundom(z):
"""
Return element of the standard fundamental domain for SL2Z
equivalent to z and a transformation in SL2Z mapping z to that
element.
EXAMPLES::
sage: from psage.ellcurve.points.chow_heegner import sl2z_rep_in_fundom
sage: z = CDF(2.17,.19)
sage: w, g = sl2z_rep_in_fundom(z)
sage: w
-2.61538461538 + 2.92307692308*I
sage: g
[ 0 1]
[-1 2]
sage: (g[0,0]*z + g[0,1])/(g[1,0]*z + g[1,1])
-2.61538461538 + 2.92307692308*I
"""
# The standard algorithm to do this is to replace z by z-n
# and z by -1/z until |Re(z)|<=1/2 and |z|>=1.
from sage.modular.all import SL2Z
gamma = SL2Z(1)
S, T = SL2Z.gens()
change = True
half = z.real().parent()(1)/2
while change:
change = False
t = z.real()
if abs(t) > half:
change = True
# |t - n| <= 1/2
# -1/2 <= t-n <= 1/2
# n - 1/2 <= t < = 1/2+n
# n <= t + 1/2 <= n + 1, so n = floor(t+1/2)
n = (t + half).floor() # avoid rounding error with 0.5
#print "t = ", t, " n = ", n
z -= n
#print "subtract", n
if hasattr(n, 'center'):
k = round(n.center())
else:
k = n
gamma *= T**k
if abs(z) < 1:
change = True
#print "invert"
z = -1/z
gamma *= S
return z, gamma**(-1)
def canonicalize_sl2z(a, g=None, eps_ratio=2):
"""
Assume that a = g(z) is in the fundamental domain.
Adjust a by applying T^(-1) or S so that a is the
canonical representative in the fundamental domain, so
a is not on the right edge, and if a is on the unit
circle, then it is on the left hand side.
Also, modify g so that the relation a = g(z) continues
to hold.
Special case: if g is None, just adjust a, ignoring g.
INPUT:
- a -- element of fundamental domain (so in ComplexField(prec))
- g -- None or element of SL2Z
- eps_ratio -- default 2; used in testing equality of the floating
point number a with 1 and 1/2. See comment in source code.
OUTPUT:
- modified a, g
"""
if g is not None:
S, T = g.parent().gens()
# There are two cases where points on boundary are identified,
# as explained in Theorem 1 of Chapter VII of Serre's "A Course
# in Arithmetic".
# Here we have to check for equality of a floating point number
# with 1/2 and with 1. In each case, exact equality of course
# almost never happens, so we must use some notion of "eps",
# and we use eps=2**(-prec//eps_ratio), where prec is the number of bits
# of precision of the parent, and eps_ratio=2 by default.
C = a.parent()
eps = 2**(-C.prec() // eps_ratio)
if abs(a.real() - C(1)/2) 0:
# points are sl2z equivalent on boundary of unit circle
a = -1/a
if g is not None: g = S*g
return a, g
def is_gamma0N_equivalent(z1, z2, N, prec):
"""
"""
w1, g1 = sl2z_rep_in_fundom(z1) # g1(z1) = w1 = canonical rep
w2, g2 = sl2z_rep_in_fundom(z2) # g2(z2) = w2 = canonical rep
#C = ComplexField(prec)
#a1 = C(w1)
#a2 = C(w2)
a1, g1 = canonicalize_sl2z(w1, g1)
a2, g2 = canonicalize_sl2z(w2, g2)
if abs(a1 - a2) > 2**(-prec):
# points are not even sl2z-equivalent, so they can't be Gamma_0(N) equivalent
return False
# TEMPORARY
g = g2**(-1) * g1
return g[1,0]%N == 0
# Now we may assume that g1(z1) = g2(z2), because of the adjustments
# made above. As a check:
assert C( (g1[0,0]*z1 + g1[0,1]) / (g1[1,0]*z1 + g1[1,1]) ) == C( (g2[0,0]*z2 + g2[0,1]) / (g2[1,0]*z2 + g2[1,1]) )
# So now we have g1(z1) = g2(z2). Let g = g2^(-1) * g1. Then g(z1) = z2.
# Suppose h is a matrix in SL2Z such that h(z1) = z2 as well.
# Then (h^(-1)*g)(z1) = h^(-1)(g(z1)) = h^(-1)(z2) = z1, so
# h^(-1)*g fixes z1. Assume that k = h^(-1)*g != 1, viewed
# as an element of PSL2Z. Then k has a fixed point in
# the upper half plane. The only elements of PSL2Z with a fixed
# point in the upper half plane are Stab(z), where
#
# * z = i, so Stab(z) generated by S
# * z = rho = exp(2*pi*i/3) so Stab(z) generated by S*T
# * z = -rhobar = exp(pi*i/3) so Stab(z) generated by T*S
#
# Assume that none of the 3 above are the case.
# Then g=h. Thus there is a matrix in Gamma_0(N) that sends
# z1 to z2 if and only if g is in Gamma_0(N), since g is the unique
# matrix in SL2Z that sends z1 to z2.
#
# In the other cases, we check the following:
# * z = i: check that neither of g, g*S in Gamma0(N)
# * z = rho: check that neither of g, g*S*T, g*(S*T)^2 in Gamma0(N)
# * z = -rhobar: check that neither of g, g*T*S, g*(T*S)^2 in Gamma0(N)
g = g2**(-1) * g1
i = C.gen()
pi = C.pi()
if g[1,0]%N == 0:
return True
S, T = g1.parent().gens()
if a1 == i:
return (g*S)[1,0]%N == 0
elif a1 == ((2*pi*i)/3).exp():
return (g*S*T)[1,0]%N == 0 or (g*S*T*S*T)[1,0]%N == 0
elif a1 == ((pi*i)/3).exp():
return (g*T*S)[1,0]%N == 0 or (g*T*S*T*S)[1,0]%N == 0
return False
class X0NPoint(object):
def __init__(self, z, N, prec):
self.z = z
self.N = N
self.prec = prec
self._C = ComplexField(prec)
def __cmp__(self, right):
assert self.N == right.N and self.prec == right.prec
if self.N == 1:
return cmp(self.sl2z_rep(), right.sl2z_rep())
if is_gamma0N_equivalent(self.z, right.z, self.N, self.prec):
return 0
return cmp(self.z, right.z)
def __repr__(self):
return "[%s]"%self._C(self.z)
#return "[%s] mod Gamma0(%s)"%(self._C(self.z), self.N)
@cached_method
def sl2z_rep(self):
a = self._C(sl2z_rep_in_fundom(self.z)[0])
b = canonicalize_sl2z(a)[0]
return b
def atkin_lehner(self, q=None):
if q is None:
# Main involution
z, N, prec = self.z, self.N, self.prec
return X0NPoint(-1/(N*z), N, prec)
else:
z, N, prec = self.z, self.N, self.prec
q = ZZ(q)
assert N%q == 0
g, x, y = xgcd(q, -N//q)
assert g == 1
# Now q*x - (N//q)*y = 1
# W_q = [q*x, y; N, q]
Wz = X0NPoint((q*x*z + y)/(N*z + q), N, prec)
#print "W_%s(%s) = %s"%(q, self, Wz)
return Wz
def __hash__(self):
return 0
# Use canonical sl2z representative.
#return hash(self.sl2z_rep())
def points_in_fiber(P, prec, low_series_prec=100, cutoff=1e-4, min_needed=None, verbose=False):
"""
INPUT:
- P -- point on elliptic curve E over a number field
- prec -- bits of precision to compute points
- index -- which low precision point to start with in upper
half plane; index=0 is the highest point, index=1 the next
highest, etc., ordered by y coordinate.
- low_series_prec -- optional low series precision used to
find low precision point in fiber.
OUTPUT:
- point z in the upper half plane to precision prec that
maps to P.
"""
E = P.curve().change_ring(QQ)
# TODO: Often elliptic logarithm totally hangs; using higher precision fixes that. :-(
# SEE trac 11767: http://trac.sagemath.org/sage_trac/ticket/11767
if verbose: print "elliptic logarithms..."
P0 = elliptic_logarithms(P, 100+prec)[0]
# Find z to low precision
f0 = phi_poly(E, low_series_prec, base_field=CDF)
if verbose: print "low precision degree = ", f0.degree()
t = walltime()
print "Finding roots in CDF of poly of degree %s"%f0.degree()
z = zeros_in_h(f0 - P0)
print "Time to find zeros of low prec series: %s seconds"%(walltime(t),)
###################
# TODO: Now refine all the zeros of f0-P0 to be full double precision
# zeros of the modular parametrization.
###################
# DO HERE!
###################
# END TODO
###################
z2 = [a for a in z if a.imag() > cutoff]
print "Will find %s points in upper half plane"%len(z2)
if min_needed:
if len(z2) < min_needed:
print "No point..."
return []
ans = []
j = 0
for tau0 in z2:
# Refine z to be a root to the requested precision
C = ComplexField(prec)
f1 = phi_poly(E, B_bound(tau0.imag(), prec=prec), base_field=C)
pi = CDF.pi(); I = CDF.gen()
q0 = CDF(exp(2*pi*I*tau0))
j += 1
if True or verbose:
print "%s/%s: refinement degree = %s (imag = %s)"%(
j, len(z2), f1.degree(), tau0.imag())
# We compute a root of f1-P0 using our assumption that the
# imaginary part is tau0.imag(). If it is ???TODO
q1 = refine_root(f1 - C(P0), C(q0), verbose=verbose)
print "q1 = ", q1
print "f1(q1)=%s, P0=%s"%(f1(q1), P0)
print "(f1-P0)(q1)=%s"%((f1-P0)(q1))
if verbose: print "refined."
tau1 = q1.log()/(2*C.pi()*C.gen())
print "first: tau0 (=%s) --> tau1 (=%s)"%(tau0, tau1)
if tau1.imag() < .8*tau0.imag():
print " **** bad. **** how good a root was tau0?"
print "(f1-C(P0))(tau0) =", (f1-C(P0))(q0)
print "(f0-P0)(tau0) =", (f0-P0)(q0)
if tau1.imag() > 0:
ans.append(tau1)
else:
print "(excluding a found point since it refined to the lower half plane)"
return ans
def index_in_mw(P):
return int((P.height()/P.curve().regulator()).sqrt().round())
def best_rational(x, contfrac=False):
if hasattr(x, 'real'):
x = x.real()
if contfrac:
return continued_fraction(x.real()).value()
else:
f = algebraic_dependency(x.real(), 1)
return f.roots(QQ, multiplicities=False)[0]
def chow_heegner_complete_fiber(E, F, x, prec,
low_series_prec=500,
prec2=53,
cutoff=1e-4,
verbose=False):
N = E.conductor()
assert F.conductor() == N
m = F.modular_degree()
P = lift_x_rational(F, x)[0]
v = points_in_fiber(P, prec=prec, low_series_prec=low_series_prec,
cutoff=cutoff, min_needed=m, verbose=verbose)
v2 = [X0NPoint(z, N=N, prec=prec2) for z in v]
v3 = [p.z for p in set(v2)]
if len(v3) > m:
raise ValueError, "found too many points (%s > %s); try increasing precision or cutoff"%(
len(v3), m)
if len(v3) < m:
raise ValueError, "not enough points (%s < %s); try increasing low_series_prec"%(
len(v3), m)
y = min(p.imag() for p in v3); y
C = ComplexField(prec)
f = phi_poly(E, B_bound(y, prec), base_field=C)
pi = C.pi(); I = C.gen()
def phi(z):
q = (2*pi*I*z).exp()
return f(q)
t1 = sum(phi(z) for z in v3)
Lambda = E.period_lattice()
t2 = Lambda.elliptic_exponential(t1)
try:
t3 = E([best_rational(t2[0]), best_rational(t2[1])])
index = int((t3.height()/E.regulator()).sqrt().round())
except Exception,msg:
print msg
t3 = None
index = None
return {'point_in_lattice':t1, 'numerical_point_on_curve':t2, 'rational_point':t3, 'index':index}
# NEXT idea: Atkin-Lehner
#
# Use Atkin-Lehner, e.g., suppose W_q(E) = eps*E for q||N.
# If eps=1, then this means W_q leaves fiber invariant, so
# we take all points we have an hit them by W_q, which is highly
# likely to double the number of points we have.
# If eps=-1, then we find points on fiber over -P (instead of P),
# and hit all those by W_q to get points on fiber ober P.
# This should massively help to make up for missing points.
def chow_heegner_atkin_lehner1(E, F, x, prec,
low_series_prec=300,
prec2=53,
cutoff=1e-4,
verbose=False):
wt = walltime(); ct = cputime()
N = E.conductor()
assert F.conductor() == N
m = F.modular_degree()
P = lift_x_rational(F, x)[0]
v = points_in_fiber(P, prec=prec, low_series_prec=low_series_prec,
cutoff=cutoff, verbose=verbose)
vplus = [X0NPoint(z, N=N, prec=prec2) for z in v]
vminus = None
v2 = copy(vplus)
v2 = list(set(v2))
if len(v2) > m:
raise ValueError, "found too many points (%s > %s); try increasing precision or cutoff"%(
len(v2), m)
for W in subsets(N.prime_divisors()):
print "W = ", W
if len(W) > 0:
if len(v2) < m:
print "not enough points (%s < %s)"%(len(v2), m)
q = prod([p**N.valuation(p) for p in W])
print "Now trying Atkin-Lehner W_%s"%q
eps = prod(F.root_number(p) for p in W)
if eps == -1 and vminus is None:
if P == -P:
vminus = vplus
else:
vminus = [X0NPoint(z,N=N,prec=prec2) for z in
points_in_fiber(-P, prec=prec, low_series_prec=low_series_prec,
cutoff=cutoff, verbose=verbose)]
v_eps = vminus
else:
v_eps = vplus
# Now for each point in v_eps, put its image under Atkin-Lehner in v2
for z in v_eps:
v2.append(z.atkin_lehner(q))
v2 = list(set(v2))
elif len(v2) > m:
raise ValueError, "found too many points (%s > %s); try increasing precision or cutoff"%(
len(v2), m)
print "Done with Atkin-Lehners: found %s of %s needed points"%(len(v2), m)
if len(v2) < m:
raise ValueError, "not enough points (%s < %s), even after using Atkin-Lehner. Trying increasing low_series_prec"%(len(v2), m)
v3 = [p.z for p in v2]
y = min(p.imag() for p in v3); y
C = ComplexField(prec)
f = phi_poly(E, B_bound(y, prec), base_field=C)
pi = C.pi(); I = C.gen()
def phi(z):
q = (2*pi*I*z).exp()
return f(q)
t1 = sum(phi(z) for z in v3)
print "Point = %s (mod lattice)"%t1
print "Computing elliptic exponential..."
Lambda = E.period_lattice()
t2 = Lambda.elliptic_exponential(t1)
try:
t3 = E([best_rational(t2[0]), best_rational(t2[1])])
index = int((t3.height()/E.regulator()).sqrt().round())
except Exception,msg:
print msg
t3 = None
index = None
return {'point_in_lattice':t1, 'numerical_point_on_curve':t2, 'rational_point':t3, 'index':index,
'walltime':walltime(wt), 'cputime':cputime(ct), 'fiber':v2, 'point_on_F':P}
def chow_heegner_atkin_lehner2(E, F, prec,
low_series_prec=300,
prec2=53,
cutoff=1e-4,
verbose=False):
wt = walltime(); ct = cputime()
N = E.conductor()
assert F.conductor() == N
m = F.modular_degree() - 1
P = F(0)
v = points_in_fiber(P, prec=prec, low_series_prec=low_series_prec,
cutoff=cutoff, verbose=verbose)
vplus = [X0NPoint(z, N=N, prec=prec2) for z in v]
v2 = copy(vplus)
v2 = list(set(v2))
if len(v2) > m:
raise ValueError, "found too many points (%s > %s); try increasing precision or cutoff"%(
len(v2), m)
for W in subsets(N.prime_divisors()):
print "W = ", W
if len(W) > 0:
if len(v2) < m:
print "not enough points (%s < %s)"%(len(v2), m)
q = prod([p**N.valuation(p) for p in W])
# Now for each point in v_eps, put its image under Atkin-Lehner in v2
for z in vplus:
v2.append(z.atkin_lehner(q))
v2 = list(set(v2))
elif len(v2) > m:
raise ValueError, "found too many points (%s > %s); try increasing precision or cutoff"%(
len(v2), m)
print "Done with Atkin-Lehners: found %s of %s needed points"%(len(v2), m)
if len(v2) < m:
raise ValueError, "not enough points (%s < %s), even after using Atkin-Lehner. Trying increasing low_series_prec"%(len(v2), m)
v3 = [p.z for p in v2]
y = min(p.imag() for p in v3); y
C = ComplexField(prec)
f = phi_poly(E, B_bound(y, prec), base_field=C)
pi = C.pi(); I = C.gen()
def phi(z):
q = (2*pi*I*z).exp()
return f(q)
t1 = sum(phi(z) for z in v3)
print "Point = %s (mod lattice)"%t1
print "Computing elliptic exponential..."
Lambda = E.period_lattice()
t2 = Lambda.elliptic_exponential(t1)
try:
t3 = E([best_rational(t2[0]), best_rational(t2[1])])
index = int((t3.height()/E.regulator()).sqrt().round())
except Exception,msg:
print msg
t3 = None
index = None
return {'point_in_lattice':t1, 'numerical_point_on_curve':t2, 'rational_point':t3, 'index':index,
'walltime':walltime(wt), 'cputime':cputime(ct), 'fiber':v2, 'point_on_F':P}
# Write Cython code to efficiently evalute a polynomial with
# ComplexField(prec) coefficients very efficiently.
#
# ANOTHER IDEA: Galois on E
#
# Another idea (not as good) -- lift Galois conjugate points from E.
# Problem with this is it requires elliptic_exponential to high precision, etc.,
# and won't scale well.
################################################################
# finding roots: not used -- use instead the dramatically faster
# chow_heegner_fast.cdf_roots_of_rdf_poly
def roots_CDF(f):
import numpy
from numpy.linalg.linalg import LinAlgError
coeffs = f.list()
numpy_array = numpy.array([complex(c) for c in reversed(coeffs)], dtype=complex)
roots = numpy.roots(numpy_array)
v = [CDF(rt) for rt in roots]
v.sort()
return v
################################################################
# Class that represents modular parametrization map. This should
# help me organize the right ideas from above, and clarify and
# fix precision issues that are killing the algorithm.
################################################################
class ModularParametrization(object):
def __init__(self, E, verbose=False):
self.verbose = verbose
self.E = E
try:
self.label = E.cremona_label()
except RuntimeError:
self.label = '%s??'%E.conductor()
def degree(self):
return self.E.modular_degree()
def cusps_over_torsion(self):
"""
Return cusps that map to a torsion point of exact order n.
"""
phi = self.E.modular_symbol_space(0).integral_period_mapping()
v = {}
for c in Gamma0(self.E.conductor()).cusps():
i = phi([c,oo])
d = i.denominator()
if v.has_key(d):
v[d].append(c)
else:
v[d] = [c]
return v
def images_of_cusps(self):
phi = self.E.modular_symbol_space(0).integral_period_mapping()
d = ZZ(1)
ims = []
for c in Gamma0(self.E.conductor()).cusps():
i = phi([c,oo])
d = d.lcm(i.denominator())
ims.append((c,i))
V = (ZZ**2)
Q = V / V.scale(d)
v = {}
for c, i in ims:
P = Q(d*i)
if v.has_key(P):
v[P].append(c)
else:
v[P] = [c]
return v
def __repr__(self):
return "Modular parametrization of %s having degree %s"%(self.label, self.degree())
def __call__(self, z):
if isinstance(z, list):
if len(z) == 0:
return []
z = [x.z for x in z]
d = max([B_bound(x.imag(), x.prec()) for x in z])
is_list = True
else:
if isinstance(z, X0NPoint):
z = z.z
z = [z]
d = B_bound(z[0].imag(), z[0].prec())
is_list = False
t = cputime()
#print "Using mod param map of degree %s to evaluate"%d
f = phi_poly(self.E, d, base_field=z[0].parent())
#print "computed phi poly"
if z[0].prec() > 53:
f = ComplexPolynomial(f)
else:
f = Polynomial_RDF_gsl(f)
#print "initialized fast poly (time: %s)"%(cputime(t))
w = []
for x in z:
q = h_to_disk(x)
tm = cputime()
w.append(f(q))
#print "time to evaluate %s"%(cputime(tm))
if not is_list:
return w[0]
return w
def _points_in_h_direct(self, z, min_imag):
prec = z.prec()
ser_prec = B_bound(min_imag, prec)
if prec < 53:
C = CDF
else:
C = z.parent()
f = phi_poly(self.E, ser_prec, base_field=C)
f -= f.parent()(z)
print f.base_ring(), C
if self.verbose: print "_points_in_h_direct -- degree: %s"%f.degree()
v = [tau for tau in zeros_in_h(f) if tau.imag() >= min_imag]
if prec < 53:
C = z.parent()
v = [C(tau) for tau in v]
return v
def _points_in_h_double(self, z, min_imag, deg1=500, max_iter=100):
z = CDF(z)
f = phi_poly(self.E, deg1, base_field=CDF) - z
v = [x for x in cdf_roots_of_rdf_poly(f) if abs(x) < 1]
print 'double: len(v) = ', len(v)
if len(v) == 0:
return []
# use actual minimum imaginary part found
w = [disk_to_h(x).imag() for x in v]
w = [x for x in w if x >= min_imag]
if len(w) == 0:
return []
min_imag = min(w)
f = phi_poly(self.E, B_bound(min_imag, 53), base_field=CDF) - z
print "deg of refinement poly = ", f.degree()
t = walltime()
g = Polynomial_RDF_gsl(f)
w = []
for b,i,err in g.newton_raphson(v, max_iter):
if abs(b) < 1 and i < max_iter:
w.append(b)
print "found %s double prec roots that refined in %s seconds"%(len(w), walltime(t))
w = throw_away_close(w,prec=45)
# Put points in upper half plane.
w = [disk_to_h(z) for z in w]
#w = [z for z in w if z.imag() >= min_imag and z.imag() <= 1000] # 1000 is effectively oo
w = [z for z in w if z.imag() >= min_imag] # 1000 is effectively oo
return w
def points_in_h(self, z, min_imag, max_iter=25, deg1=500, max_iter1=1000, equiv_prec=None):
"""
Return as points tau in open upper half plane with tau.imag()
>= min_imag to precision z.prec() that all map to z via the
modular parametrization map.
Doesn't use Atkin-Lehner's.
INPUT:
- z -- complex number, view as a point in C / Lambda, i.e.,
the elliptic logarithm of a point on the curve self.E
- min_imag -- positive real number
- deg1 -- degree used for first double precision root
finding; if too small then some points will be missed.
- max_iter -- maximum number of iterations used in
newton-raphson calls
OUTPUT:
- list of complex numbers
"""
C = z.parent()
N = self.E.conductor()
if C != CDF and not is_ComplexField(C):
raise ValueError, "z must be in complex field"
v = self._points_in_h_double(z, min_imag=min_imag, deg1=deg1, max_iter=max_iter1)
prec = z.prec()
if equiv_prec is None:
equiv_prec = prec//4
if prec <= 53:
return [X0NPoint(C(x),N,prec=equiv_prec) for x in v]
# refine and keep good roots
if len(v) == 0:
return []
print "len(v) =", len(v)
t = walltime()
f = phi_poly(self.E, B_bound(min_imag, prec), base_field=C) - z
w = []
for b,i,err in newton_raphson(f, [h_to_disk(C(a)) for a in v], max_iter, max_err=C(2)**(2-prec)):
#print 'b2', b,i,err
if abs(b) < 1 and i < max_iter:
w.append(b)
w = throw_away_close(w,prec=prec-10)
w = [disk_to_h(z) for z in w]
w = [z for z in w if z.imag() >= min_imag and z.imag() <= 1] # above 1 --> is point at oo
w = [X0NPoint(z, N, prec=prec//2) for z in w]
w = list(set(w))
w.sort()
print "found %s roots to prec %s that refined in %s seconds"%(len(w), prec, walltime(t))
return w
def disk_to_h(q):
K = q.parent()
return q.log() / (2*K.pi()*K.gen())
def h_to_disk(z):
K = z.parent()
return (2*K.pi()*K.gen()*z).exp()
def throw_away_close(v, prec):
"""
Return list of 'distinct' elements of v, where we consider two
elements close if they agree when coerced down to prec bits
of precision.
"""
C = ComplexField(prec)
class W:
def __init__(self, x):
self.x = x
self.y = C(x)
def __hash__(self):
return hash(self.y)
def __cmp__(self, right):
return cmp(self.y, right.y)
w = [a.x for a in set([W(x) for x in v])]
w.sort()
return w
def newton_raphson(f, x, max_iter=1000, max_err=1e-14):
"""
Use the given number of steps of the Newton-Raphson algorithm
to refine the given input roots of f.
If more than one root as input, returns list of refined roots
and number of iterations and error.
"""
C = f.base_ring()
if C == CDF:
return Polynomial_RDF_gsl(f).newton_raphson(x, max_iter=max_iter, max_err=max_err)
if C.prec() < 53:
t = Polynomial_RDF_gsl(f).newton_raphson(x, max_iter=max_iter, max_err=max_err)
return (C(t[0]), t[1], t[2])
if not isinstance(x, list):
x = [x]
f_prime = f.derivative()
# much faster evaluation versions
f = ComplexPolynomial(f)
f_prime = ComplexPolynomial(f_prime)
ans = []
for root in x:
last_root = C(root)
for i in range(max_iter):
root = root - f(root)/f_prime(root)
err = abs(last_root - root)
if err <= max_err:
break
last_root = root
ans.append((root, i+1, err))
return ans
def chow_heegner_atkin_lehner3(E, F, x, prec,
low_series_prec=300,
prec2=53,
cutoff=1e-4,
verbose=False):
wt = walltime(); ct = cputime()
N = E.conductor()
assert F.conductor() == N
m = F.modular_degree()
P = lift_x_rational(F, x)[0]
v = points_in_fiber(P, prec=prec, low_series_prec=low_series_prec,
cutoff=cutoff, verbose=verbose)
vplus = [X0NPoint(z, N=N, prec=prec2) for z in v]
vminus = None
v2 = copy(vplus)
v2 = list(set(v2))
if len(v2) > m:
raise ValueError, "found too many points (%s > %s); try increasing precision or cutoff"%(
len(v2), m)
for W in subsets(N.prime_divisors()):
print "W = ", W
if len(W) > 0:
if len(v2) < m:
print "not enough points (%s < %s)"%(len(v2), m)
q = prod([p**N.valuation(p) for p in W])
print "Now trying Atkin-Lehner W_%s"%q
eps = prod(F.root_number(p) for p in W)
if eps == -1 and vminus is None:
if P == -P:
vminus = vplus
else:
vminus = [X0NPoint(z,N=N,prec=prec2) for z in
points_in_fiber(-P, prec=prec, low_series_prec=low_series_prec,
cutoff=cutoff, verbose=verbose)]
v_eps = vminus
else:
v_eps = vplus
# Now for each point in v_eps, put its image under Atkin-Lehner in v2
for z in v_eps:
v2.append(z.atkin_lehner(q))
v2 = list(set(v2))
elif len(v2) > m:
raise ValueError, "found too many points (%s > %s); try increasing precision or cutoff"%(
len(v2), m)
print "Done with Atkin-Lehners: found %s of %s needed points"%(len(v2), m)
if len(v2) < m:
raise ValueError, "not enough points (%s < %s), even after using Atkin-Lehner. Trying increasing low_series_prec"%(len(v2), m)
v3 = [p.z for p in v2]
y = min(p.imag() for p in v3); y
C = ComplexField(prec)
f = phi_poly(E, B_bound(y, prec), base_field=C)
pi = C.pi(); I = C.gen()
def phi(z):
q = (2*pi*I*z).exp()
return f(q)
t1 = sum(phi(z) for z in v3)
print "Point = %s (mod lattice)"%t1
print "Computing elliptic exponential..."
Lambda = E.period_lattice()
t2 = Lambda.elliptic_exponential(t1)
try:
t3 = E([best_rational(t2[0]), best_rational(t2[1])])
index = int((t3.height()/E.regulator()).sqrt().round())
except Exception,msg:
print msg
t3 = None
index = None
return {'point_in_lattice':t1, 'numerical_point_on_curve':t2, 'rational_point':t3, 'index':index,
'walltime':walltime(wt), 'cputime':cputime(ct), 'fiber':v2, 'point_on_F':P}
class NumericalPoint(object):
def __init__(self, P, eps, tau=None):
self.P = P
self.eps = eps
self.tau = tau
def curve(self):
return self.P.curve()
def rational_curve(self):
return EllipticCurve(QQ,[int(a.real()) for a in self.P.curve().a_invariants()])
def index(self, **kwds):
"""
Index of point in the Mordell-Weil group modulo torsion.
"""
E = self.rational_curve()
assert E.rank() == 1
h = self.best_approx(**kwds).height()
return int((h / E.regulator()).sqrt().round())
def best_approx(self, known_bits=None, infinity=1e15):
if known_bits is None:
known_bits = int(self.P[0].prec()*.7) # heuristic
E = self.rational_curve()
if abs(self.P[0]) >= infinity or abs(self.P[1]) >= infinity or abs(self.P[2]) < 1.0/infinity:
return E(0)
x = algebraic_dependency(self.P[0], 1, known_bits=known_bits).roots(QQ, multiplicities=False)[0]
v = E.lift_x(x, all=True)
if len(v) == 0:
raise ValueError, "no rational point on curve that approximates this point"
m = [(abs(w[1] - self.P[1]), w) for w in v]
m.sort()
return m[0][1]
def close_points(self, search_bound, eps=1e-3):
E = self.rational_curve()
g = E.gens()
if len(g) == 0:
P = E(0)
search_bound=0
else:
P = g[0]
T = E.torsion_points()
v = []
for n in range(-search_bound, search_bound+1):
Q = n*P
for t in T:
R = Q + t
if abs(R[0] - self[0]) < eps and abs(R[1] - self[1]) < eps:
v.append(R)
return v
def __getitem__(self, i):
return self.P.__getitem__(i)
def __hash__(self):
return int(0)
def __add__(self, right):
if isinstance(right, int) and right == 0:
return self
return NumericalPoint(self.P + self.P.curve().point((right.P[0],right.P[1],right.P[2]),check=False), self.eps)
def __radd__(self, left):
if isinstance(left, int) and left == 0:
return self
raise NotImplementedError
def __rmul__(self, left):
return NumericalPoint(left*self.P, self.eps)
def __mul__(self, right):
return NumericalPoint(self.P*right, self.eps)
def __cmp__(self, right):
if max([abs(self.P[i]-right.P[i]) for i in range(3)]) < min(self.eps, right.eps):
return 0
return cmp(self.P, right.P)
def __repr__(self):
return repr(self.P)
def identify(P, search=200, eps=1e-5, infinity=1e12):
if max(abs(P[0]),abs(P[1])) >= infinity or P.P == 0:
return P.rational_curve()(0), 0
m = P.close_points(search, eps)
if len(m) != 1:
raise ValueError, "unable to identify rational point"
else:
P = m[0]
return P, index_in_mw(P)
class ChowHeegner2(object):
"""
Low precision strategy for finding Chow-Heegner points.
"""
def __init__(self, E, F):
if isinstance(E, str): E = EllipticCurve(E)
if isinstance(F, str): F = EllipticCurve(F)
if E.conductor() != F.conductor():
print "*"*70
print "BIG FAT WARNING: Conductors are different -- probably get total nonsense!"
print "*"*70
self.E = E
self.F = F
self.fE = ModularParametrization(E)
self.fF = ModularParametrization(F)
def __repr__(self):
return "Chow-Heegner(%s, %s)"%(label(self.E), label(self.F))
def points_over_F(self, z=CDF(0.1), deg1=500, min_imag=1e-3, max_iter1=100, equiv_prec=None,
use_atkin_lehner=False):
z = complexify(z)
assert z.prec() >= 53
v = list(set(self.fF.points_in_h(z, deg1=deg1, min_imag=min_imag, max_iter1=max_iter1, equiv_prec=equiv_prec)))
if not use_atkin_lehner:
return v
F = self.F
N = F.conductor()
for i, p in enumerate(N.prime_divisors()):
if F.root_number(p) == 1:
q = p**N.valuation(p)
v = list(set(v + [x.atkin_lehner(q) for x in v]))
v.sort()
return v
def _points_on_E(self, z=CDF(0.1), deg1=500, min_imag=1e-3, eps=1e-3, max_iter1=100):
z = complexify(z)
assert z.prec() >= 53
v = self.points_over_F(z=z, deg1=deg1, min_imag=min_imag, max_iter1=max_iter1)
w = []
if z == 0:
w.append(NumericalPoint(self.E(0),eps=eps))
Lambda = self.E.period_lattice()
for tau in v:
P = Lambda.elliptic_exponential(self.fE(tau))
w.append(NumericalPoint(P, eps=eps, tau=tau))
w = list(set(w))
w.sort()
return w
def _fiber(self, z, target, min_equiv_prec=10,
use_atkin_lehner=False, # horrible
max_prec=1e10, deg1=200, min_imag=1e-3, equiv_prec=10):
if equiv_prec target:
print "too many points: changing params"
if equiv_prec + 1 <= z.prec()//2:
equiv_prec += 1
recompute = False
else:
if z.prec()+5 <= max_prec:
z = ComplexField(z.prec()+5)(z)
elif len(w) < target:
print "not enough points: changing params"
# increase parameters
if equiv_prec - 1 >= min_equiv_prec:
equiv_prec -= 1
recompute = False
else:
deg1 += 300
if min_imag >= 2*(1e-6):
min_imag /= 2
if z.prec()+5 <= max_prec:
z = ComplexField(z.prec()+5)(z)
elif len(w) == target:
print "*** W00t! found all needed points ***"
return w, deg1, min_imag, equiv_prec, z.prec()
def point_on_E_fiber_cardinality(self, z=CDF(.1), min_equiv_prec=10,
use_atkin_lehner=False, # horrible
max_prec=1e10):
tt = cputime()
mF = self.fF.degree()
print "modular degree of F (target fiber size) =", mF
target = mF
w, deg1, min_imag, equiv_prec, prec = self._fiber(
z, target, min_equiv_prec, use_atkin_lehner=use_atkin_lehner,max_prec=max_prec)
print "mapping them to E..."
t = cputime()
m = self.fE(w)
print "mapped to E in %s seconds"%cputime(t)
P = NumericalPoint(self.E.elliptic_exponential(sum(m)), 1e-4)
return {'P':P, 'prec':prec, 'deg1':deg1, 'min_imag':min_imag,
'equiv_prec':equiv_prec, 'cputime':cputime(tt)}
def point_on_E_target_cardinality(self, z=CDF(0.1)):
tt = cputime()
deg1=300
min_imag=1e-3
mF = self.fF.degree()
print "modular degree of F =", mF
# assume fiber is generic and z doesn't correspond to a 2 torsion point.
n = len([p for p in self.E.conductor().prime_divisors() if
self.E.root_number(p) == 1 and self.F.root_number(p) == 1])
if target is None and z != 0:
n = self._atkin_lehner_factor()
if mF % n == 0:
target = mF // n
else:
target = mF
print "WARNING: mF (=%s) not divisible by (#agreeing w)=%s"%(mF, n)
multiplier = n
else:
multiplier = 1
print "target number of points on E =", target
w = []
while len(w) < target:
b = len(w)
print "** deg1 = %s, min_imag = %s, len(w) = %s, target = %s"%(deg1, min_imag, len(w), target)
tm = cputime()
for t in self._points_on_E(z, deg1=deg1, min_imag=min_imag):
w.append(t)
w = list(set(w))
print "** found %s new points in %s seconds"%(len(w)-b, cputime(tm))
if len(w) < target:
# increase parameters
deg1 += 300
min_imag /= 2
elif len(w) == target:
print "*** W00t! found all needed points ***"
else:
raise RuntimeError, "found too many points (%s > %s) -- precision bug"%(len(w),target)
P = multiplier * sum(w)
return {'P':P, 'w':w, 'cputime':cputime(tt),
'deg1':deg1, 'min_imag':min_imag, 'z':z}
def point_on_E_using_0_fiber(self, prec=53, min_equiv_prec=10, max_prec=1e10, **kwds):
tt = cputime()
# 1. Find the cusps of X_0(N) that map to 0 in F, using modular symbols.
C = self.fF.cusps_over_torsion()[1]
print "%s cusps map to 0: %s"%(len(C), C)
# 2. Find mod_deg(F) - len(C) more distinct affine points on X_0(N) numerically.
target = self.fF.degree() - len(C)
w, deg1, min_imag, equiv_prec, prec = self._fiber(ComplexField(prec)(0), target, min_equiv_prec, max_prec=max_prec, **kwds)
# 3. Map down non-cuspidal points and sum
m = self.fE(w)
P = NumericalPoint(self.E.elliptic_exponential(sum(m)), 1e-4)
return {'P':P, 'prec':prec, 'deg1':deg1, 'min_imag':min_imag,
'equiv_prec':equiv_prec, 'cputime':cputime(tt)}
def _atkin_lehner_factor(self):
# We assume fiber is generic and z doesn't correspond to element of F[2]
# Let G be the subgroup of the Atkin-Lehner group W consisting of all
# elements of W that fix both E and F. Then I think we can divide
# the target cardinality by 2^#G.
N = self.E.conductor()
wE, wF = [[X.root_number(p) for p in N.prime_divisors()] for X in [self.E, self.F]]
n = 1
for s in subsets(range(len(wE))):
if len(s) > 0:
#print prod(wE[i] for i in s), prod(wF[i] for i in s)
if prod(wE[i] for i in s) == 1 and prod(wF[i] for i in s) == 1:
n += 1
return n
def point_on_E_target_cardinality2(self, z=CDF(0.1), deg1=500, min_imag=1e-4,
target=None, multiplier=1,
number_to_stabilize=5, w=None):
tt = cputime()
mF = self.fF.degree()
print "modular degree of F =", mF
if target is None and z != 0:
n = self._atkin_lehner_factor()
if mF % n == 0:
target = mF // n
else:
target = mF
print "WARNING: mF (=%s) not divisible by (#agreeing w)=%s"%(mF, n)
multiplier = n
print "target number of points on E =", target
base_points = [z]
Omega = self.F.period_lattice().basis(z.prec())[0]
if w is None:
w = []
else:
w = list(w)
n = 0
same_in_a_row = 0
stabilize = False
while len(w) < target:
b = len(w)
print "** deg1 = %s, min_imag = %s, len(w) = %s, target = %s"%(deg1, min_imag, len(w), target)
tm = cputime()
for t in self._points_on_E(z, deg1=deg1, min_imag=min_imag):
w.append(t)
w = list(set(w))
print "** found %s new points in %s seconds"%(len(w)-b, cputime(tm))
if len(w) == b:
same_in_a_row += 1
if same_in_a_row >= number_to_stabilize:
print "** same result %s times in a row, so 'stabilized'"%same_in_a_row
stabilize = True
break
else:
same_in_a_row = 0
if len(w) < target:
# changing base point
n += 1
z += ((-1)**(n+1)) * n*Omega
base_points.append(z)
print "z (=%s) |--> %s"%(base_points[-1], z)
elif len(w) == target:
print "*** W00t! found all needed points ***"
else:
raise RuntimeError, "found too many points (%s > %s) -- precision bug"%(len(w),target)
P = multiplier * sum(w)
return {'P':P, 'w':w, 'cputime':cputime(tt),
'deg1':deg1, 'min_imag':min_imag, 'base_points':base_points,
'E':label(self.E), 'F':label(self.F),
'stabilize':stabilize,
'multiplier':multiplier}
def point_on_E_fiber_cardinality2(self, z=CDF(.1),
equiv_prec=10, deg1=500, min_imag=1e-4,
fiber=None,
map_to_E=True,
number_to_stabilize=6):
z = complexify(z)
tt = cputime()
if fiber is None:
fiber = []
else:
fiber = list(fiber)
base_points = [z]
mF = self.fF.degree()
print "modular degree of F (target fiber size) =", mF
target = mF
Omega = self.F.period_lattice().basis(z.prec())[0]
n = 0
if z == 0:
# reduce target by number of cusps that map to 0
C = self.fF.cusps_over_torsion()[1]
target -= len(C)
print "reducing target by %s cusps to %s (TODO *** result may be off by torsion! *** )"%(len(C), target)
while len(fiber) < target:
b = len(fiber)
print "** deg1 = %s, min_imag = %s, equiv_prec = %s, len(fiber) = %s, target = %s"%(deg1, min_imag, equiv_prec, b, target)
tm = cputime()
for t in self.points_over_F(z, min_imag=min_imag, deg1=deg1, equiv_prec=equiv_prec):
fiber.append(t)
fiber = list(set(fiber))
print "** found %s new points in %s seconds"%(len(fiber)-b, cputime(tm))
if len(fiber) == b:
same_in_a_row += 1
if same_in_a_row >= number_to_stabilize:
print "*"*80
print "** ERROR: same result %s times in a row, so 'stabilized' (GIVING UP!)"%same_in_a_row
print "** Result surely nonsense."
print "*"*80
stabilize = True
break
else:
same_in_a_row = 0
if len(fiber) < target:
print '*'*80
# change base point
n += 1
z += ((-1)**(n+1)) * n*Omega
base_points.append(z)
print "z (=%s) |--> %s"%(base_points[-1], z)
elif len(fiber) == target:
print '+'*80
print "*** W00t! found all needed points ***"
else:
print '-'*80
raise RuntimeError, "found too many points (%s > %s) -- try (good) increasing precision of z (now=%s) or (bad) *decreasing* equiv_prec (now=%s)"%(
len(fiber),target,z.prec(),equiv_prec)
data = {'base_points':base_points, 'deg1':deg1, 'min_imag':min_imag,
'equiv_prec':equiv_prec, 'cputime':cputime(tt),
'fiber':fiber}
if map_to_E:
print "mapping them to E..."
t = cputime()
m = self.fE(fiber)
print "mapped to E in %s seconds"%cputime(t)
s = sum(m)
if isinstance(s, int) and s == 0:
P = NumericalPoint(self.E(0), 1e-4)
else:
P = NumericalPoint(self.E.elliptic_exponential(s), 1e-4)
data['P'] = P
return data
def point_on_E_2(self,
z=CDF(.1), equiv_prec=10, deg1=500, min_imag=1e-4,
number_to_stabilize=6,
point_eps=1e-3,
fiber=None, points_on_E=None):
"""
Combine both fiber and card image methods.
1. Using the given parameters, compute points in fiber in h.
2. If enough points in h, done.
3. If not, map those points down to E.
4. If enough points in E, done.
5. Go to 1, adding to what we have so far.
"""
input_params = locals()
warnings = []
z = complexify(z)
tt = cputime()
if fiber is None:
fiber = []
else:
fiber = list(fiber)
if points_on_E is None:
points_on_E = []
else:
points_on_E = list(points_on_E)
base_points = [z]
mF = self.fF.degree()
print "modular degree of F (target fiber size) =", mF
target_fiber = mF
Omega = self.F.period_lattice().basis(z.prec())[0]
n = 0
if z == 0:
# reduce target by number of cusps that map to 0
C = self.fF.cusps_over_torsion()[1]
if len(C) > 0:
target_fiber -= len(C)
s = "reducing target by %s cusps to %s (TODO *** result may be off by torsion! *** )"%(len(C), target_fiber)
print s
warnings.append(s)
n = self._atkin_lehner_factor()
if mF % n == 0:
target_E = mF // n
else:
target_E = mF
s = "WARNING: mF (=%s) not divisible by (#agreeing w)=%s"%(mF, n)
print s
warnings.append(s)
multiplier_on_E = n
reports = []
timestamp = cputime()
def report():
s = "** %.1fs: z=%s, deg1=%s, min_imag=%e, equiv_prec=%s, len(fiber)=%s<=%s, len(points_on_E)=%s<=%s"%(
cputime(timestamp), z, deg1, min_imag, equiv_prec, len(fiber), target_fiber, len(points_on_E), target_E)
reports.append(s)
return s
fail = False
while len(fiber) < target_fiber and len(points_on_E) < target_E:
b = len(fiber)
tm = cputime()
print report()
new_points = self.points_over_F(z, min_imag=min_imag, deg1=deg1, equiv_prec=equiv_prec)
last_fiber = set(fiber)
for t in new_points:
fiber.append(t)
fiber = list(set(fiber))
print "** found %s new points in fiber in %s seconds"%(len(fiber)-b, cputime(tm))
Lambda = self.E.period_lattice()
if len(fiber) == b:
same_in_a_row += 1
if same_in_a_row >= number_to_stabilize:
print "*"*80
print fail
print "*"*80
fail = ('stabilize', "** ERROR: same result %s times in a row, so GIVING UP!"%same_in_a_row)
break
else:
same_in_a_row = 0
# Now add points to curve.
b0 = len(points_on_E)
for tau in set(fiber).difference(last_fiber):
# map it down.
P = Lambda.elliptic_exponential(self.fE(tau))
points_on_E.append(NumericalPoint(P, eps=point_eps, tau=tau))
points_on_E = list(set(points_on_E))
print "** gave rise to %s new points on E in %s seconds"%(len(points_on_E)-b0, cputime(tm))
if len(fiber) < target_fiber and len(points_on_E) < target_E:
print '*'*80
# change base point
n += 1
z += ((-1)**(n+1)) * n*Omega
base_points.append(z)
elif len(fiber) == target_fiber:
print '+'*80
print "*"*70
print "** W00t! found enough points on fiber! ***"
print "*"*70
elif len(points_on_E)==target_E:
print '+'*80
print "*"*70
print "** W00t! found enough points on E! ***"
print "*"*70
elif len(fiber) > target_fiber:
print '-'*80
fail = ('fiber', "found too many points (%s > %s) in fiber -- try (good) increasing precision of z (now=%s) or (bad) *decreasing* equiv_prec (now=%s)"%(
len(fiber),target_fiber,z.prec(),equiv_prec))
print fail
break
elif len(points_on_E) > target_E:
print '-'*80
fail = ('points_on_E', "found too many points (%s > %s) on E -- try (good) increasing precision of z (now=%s) or (bad) *decreasing* equiv_prec (now=%s)"%(
len(points_on_E), target_E, z.prec(), equiv_prec))
print fail
break
print report()
data = {'base_points':base_points, 'deg1':deg1, 'min_imag':min_imag,
'equiv_prec':equiv_prec, 'cputime':cputime(tt),
'fiber':fiber,
'points_on_E':points_on_E,
'number_to_stabilize':number_to_stabilize,
'reports':reports,
'warnings':warnings,
'input_params':input_params
}
if not fail:
if len(fiber) == target_fiber:
print "mapping them to E..."
t = cputime()
m = self.fE(fiber)
print "mapped to E in %s seconds"%cputime(t)
s = sum(m)
if isinstance(s, int) and s == 0:
P = NumericalPoint(self.E(0), 1e-4)
else:
P = NumericalPoint(self.E.elliptic_exponential(s), 1e-4)
data['P'] = P
data['use'] = 'fiber'
elif len(points_on_E) == target_E:
P = multiplier_on_E * sum(points_on_E)
data['P'] = P
data['multiplier_on_E'] = multiplier_on_E
data['use'] = 'points_on_E'
else:
assert False, "bug"
try:
ident = identify(P)
except ValueError:
ident = ("failed to identify", "?")
data['id'] = ident
else:
data['use'] = 'fail'
data['fail'] = fail
return data
############################################
def complexify(z):
C = parent(z)
if C == CDF or is_ComplexField(C):
return z
try:
prec = z.prec()
except AttributeError:
prec = 53
C = CDF if prec==53 else ComplexField(prec)
return C(z)
##############################
# Another idea:
# Now that we can correctly divide up points, and figure out which cusps map to 0,
# implement the "obvious function" in the case when 0 is unramified, i.e.,
# (1) find all cusps that map to 0
# (2) find exactly the right number of remaining points in upper half plane.
# (3) map only the remaining points down.
# This might work in some cases...?
######
##################################################################
# TABLES
##################################################################
def rank1_curves(levels):
for E in cremona_optimal_curves(levels):
if E.rank() == 1:
yield E
@parallel
def table_row(E, F, compute=True):
Elabel = E.cremona_label()
Flabel = F.cremona_label()
row = [Elabel, {'F':Flabel, 'degE':E.modular_degree(), 'degF':F.modular_degree()}]
print row
if not compute:
return row
filename = 'data/point_on_E_2/%s-%s.sobj'%(Elabel, Flabel)
if not os.path.exists(os.path.split(filename)[0]):
os.makedirs(os.path.split(filename)[0])
if os.path.exists(filename):
data = load(filename)
print data
else:
ch = ChowHeegner2(E, F)
print "*"*70
print ch
data = ch.point_on_E_2()
print data
save(data, filename)
def table(levels, compute=True):
for E in rank1_curves(levels):
for row in table_curve(E, compute=compute):
yield row
def table_curves(E):
for F in cremona_optimal_curves([E.conductor()]):
if F != E:
yield (E, F)
def table(levels, compute=True, parallel=True):
v = sum([list(table_curves(E)) for E in rank1_curves(levels)], [])
if not compute:
for X in v:
yield X
else:
if parallel:
for X in table_row(v):
print X
else:
for E,F in v:
print table_row(E,F,parallel=False)
## print table
def plain_table(outfile = 'table.txt'):
"""
Iterate over the pickle files in the current directory,
and output a simple easy-to-read file summarizing the
contents. The columns are:
E F index degE degF rankF cputime
"""
out = open(outfile,'w')
out.write('E\tF\tindex\tdegE\tdegF\trankF\tcputime\n')
for v in os.listdir(os.curdir):
if not v.endswith('sobj'): continue
X = load(v)
w = v.split('-')
Elabel = w[0]
Flabel = w[1].split('.')[0]
E = EllipticCurve(Elabel)
F = EllipticCurve(Flabel)
if X['use'] == 'fail':
ind = '?'
else:
ind = X['id'][1]
tm = X['cputime']
s = '%s\t%s\t%s\t%s\t%s\t%s\t%.1f'%(
Elabel[:-1], Flabel[:-1], ind, E.modular_degree(),
F.modular_degree(), F.rank(), tm)
print s
out.write(s+'\n')