""""\Teichmuller curves in genus three and \\ just likely intersections in $\IG_m^n\times\IG_a^n$."
by Matt Bainbridge, Philipp Habegger, Martin M\"oller
Description: Sage code to varify computational claims of the Proof of Proposition 9.6. Any failure leads to an assertion falt.
Usage: Create an empty directory called "Data" to store intermediate calculations. This directory
should be located in the "working_directory" whose path is specified below (delete the contents of this directory if you want to start the computations from scratch). Then run "sage g3fin_ch9.sage".
Works with sage 6.6.
"""
start_time = walltime()
#set to directory containing the program files
working_directory = './'
file = working_directory + 'Data/'
from os import mkdir
"""
make data directory if it doesn't exist.
To restart the computation from scratch, delete everything in this directory.
"""
try:
mkdir(file)
except:
pass
"""
load(working_directory + 'subspaces.sage')
"""
def build_all_subspaces(matrix_list, eqlist):
"""
build_all_subspaces builds the list of matrices in the first part of the algorithm of Section 5
matrix_list is a list of initial subspaces of the ambient Z^N. We are looking for tori parallel
to one whose associated subspace is contained in one in this list.
eqlist is the list of equations cutting out the variety we are searching for subtori of. This
module only uses the list of exponents of these equations.
the output subspaces is a list of lists of matrices. Subspaces[i] are the subspaces which are
codimension i in some initial matrix in matrix_list.
"""
R = matrix_list[0].rank() #all initial matrices should have the same rank
assert([M1.rank() == R for M1 in matrix_list])
subspace_list=[matrix_list]
for i in range(R-1):
subspace_set = []
for N in subspace_list[-1]:
new_spaces = build_subspaces(N, eqlist)
for space in new_spaces:
if space not in subspace_set:
subspace_set.append(space)
subspace_list.append(subspace_set)
return subspace_list
def build_subspaces(M,eqlist):
"""
M is an integer matrix whose rows span a subspace of Z^N.
eqlist is a list of polynomials.
build_subspaces outputs the potential codimension-one subspaces of span(M) which are candidates
for a torus in the variety cut out by eqlist.
"""
subspaces = []
for eq in eqlist:
projection_list = [orthogonal_projection(vector(e), M) for e in eq.exponents()] #list of projections of exponent vectors to M
L = len(projection_list)
for i in range(0, L):
for j in range(i+1, L):
v = projection_list[i] - projection_list[j]
if v != 0:
subspace = orthogonal_complement(v, M)
subspace.set_immutable()
if subspace not in subspaces:
subspaces.append(subspace)
return subspaces
def orthogonal_complement(v, M):
"""
returns basis of the orthogonal complement of v in the row-span of M
Assumes: v is contained in the row-span and v!0 =.
"""
assert v != 0
N = matrix([v] + M.rows())
assert N.rank() == M.rank()
N = N.gram_schmidt()[0]
assert N.rows()[0] == v #gram-schmidt shouldn't change the first row.
return make_integral(N[1:].echelon_form()) #return the rest of the rows in echelon form. Since the first row is v, and the rows are orthogonal, these rows span the orthogonal complement of v in span(M)
def orthogonal_projection(w, M):
"""
orthogonal projection of w onto the row span of M
"""
M = M.gram_schmidt()[0] #orthogonal basis of the span of M
return sum([(v.dot_product(w) / v.dot_product(v) ) * v for v in M.rows() ] )
def make_integral(N):
"""
N is a rational matrix. Returns integer matrix obtained by multiplying each row of N by a constant so that the resulting row is integral with no common divisor.
"""
N1 = matrix([ v*v.denominator() for v in N.rows()]).change_ring(ZZ)
N2 = matrix([v/gcd(v) for v in N1.rows()]).change_ring(ZZ)
return N2
def get_equation_list(eq, M):
"""
eq is a polynomial, M an integral matrix.
returns the list of equations expressing that a torus-translate parallel to the torus defined by
M is contained in the hypersurface cut out by eq. This is the system of polynomials p_J defined
in the last displayed equation of Section 5.
"""
eqdict = eq.dict()
d = {}
for e in eqdict.keys():
ve = vector(e)
w = orthogonal_projection(ve, M)
wt = tuple(w.list())
#print e, w
d[wt] = d.get(wt, 0) + eqdict[e]*prod([a[i]^e[i] for i in range(0, 9)])
#assert(sum(d.values()) == eq)
#print d
return [x for x in d.values()]
"""
def get_all_equations_old(eqlist, M):
equations = []
for eq in eqlist:
equations.extend(get_equation_list(eq, M))
return Pa.ideal(equations)
"""
matrix_list = []
"""
Enter here the list of matrices determining the initial tori. The algorithm searches for
torus-translates parallel to subtori of one in this list.
"""
"""
This is M3 in the paper
"""
M = matrix([[-1,0,0,1,0,0,0,0,0],[0,-1,0,0,0,0,0,1,0],[0,0,0,0,0,1,0,0,-1]])
M.set_immutable()
matrix_list.append(M)
"""
This is M1 in the paper
"""
M = matrix([[1,1,1,1,1,1,0,0,0],[0,1,1,0,1,1,0,1,1],[0,0,1,0,0,1,0,0,0]])
M.set_immutable()
matrix_list.append(M)
"""
This is M2 in the paper
"""
M = matrix([[1,1,1,1,1,1,0,0,0],[0,1,1,0,1,1,0,1,1],[0,0,0,0,0,1,0,0,1]])
M.set_immutable()
matrix_list.append(M)
"""
The variables pi and qi correspond to xi and yi in the paper. The indicies have been shifted by 1
to accomodate the usual python convention of arrays starting with 0. y3 corresponds to the "free"
zero z4 (the other three are 0, 1, infinity.
The rij correspond to the cross-ratios R{i+1}{j+1}
"""
P. = PolynomialRing(QQ)
Q. = PolynomialRing(QQ)
y = [Infinity,0,1,y3]
p = [p0, p1, p2]
q = [q0, q1, q2]
r = [[],[r10,r11,r12],[r20,r21,r22],[r30,r31,r32]]
def CR(a, b, c, d): #cross ratio [a,b,c,d]
if a == Infinity:
return (b-d)/(b-c)
else:
return (a-c)*(b-d)/((a-d)*(b-c))
def R(i,j,k): #corresponds to R_ijk from the paper
return CR(y[i], y[j], p[k], q[k])
"""
f represents the cross-ratio morphism from M_0,10 to the subvariety Y of A_red
"""
f = FractionField(Q).hom([R(0,1,0), R(0,2,0), R(0,3,0), R(0,1,1), R(0,2,1), R(0,3,1), R(0,1,2), R(0,2,2), R(0,3,2)], FractionField(P))
"""
f is birational. Here is its inverse:
"""
g = FractionField(P).hom([((r20-1)/(r20-r10)*r10 - r30*(r20-1)/(r20-r10))/(1-r30),(r20-1)/(r20-r10),(r21-1)/(r21-r11),(r22-1)/(r22-r12), (r20-1)/(r20-r10)*r10, (r21-1)/(r21-r11)*r11, (r22-1)/(r22-r12)*r12], FractionField(Q))
"""
Note that g does not involve r31 or r32. This means that we can regard g as either induced by a
rational map (Y \subset CR^\red)--> M_{0,10}, or a rational map CR^\min --> M_{0,10}. These maps
are related by the projection p_\min: CR^\red --> CR^\min.
We check below that g is birational
"""
"""
Yideal is the ideal I from the paper, cutting out the closure of \cal{Y}
"""
try:
Yideal = load(file + 'Yideal')
except IOError:
eq0 = (CR(1, r10, r20, r30) - CR(1, r11, r21, r31)).numerator() #Eq(4,1,2)
eq1 = (CR(1, r10, r20, r30) - CR(1, r12, r22, r32)).numerator() #Eq(4,1,3)
eq2 = (CR(1, r11, r21, r31) - CR(1, r12, r22, r32)).numerator() #Eq(4,2,3)
Yideal = Q.ideal([eq0, eq1,eq2])
save(Yideal, file + 'Yideal')
assert(Yideal.is_prime()) #check that Y is Q-irreducible
assert(f(Yideal) == FractionField(P).ideal(0)) #check that the image of f is contained in Y
Yring = Q.quotient(Yideal) #the coordinate ring of Y
"""
Lets verify that f is birational onto Y with inverse g.
"""
FPid = FractionField(P).hom(P.gens(), FractionField(P)) #identity homomorphism of FractionField(P)
pi = Yring.cover() #caononical morphism Q-->Yring
pifrac = FractionField(Q).hom(Yring.gens(), FractionField(Yring)) #induced morphism on fraction fields
assert(f*g == FPid) #f*g is the identity on the function field of M_{0,10}
for i in range(1,3): #check that g*f is the identity on the coordinates r10, r11, r12, r20, r21, r22, r30 of A^\min
for j in range(0,3):
assert((g*f)(r[i][j]) == r[i][j])
assert((g*f)(r[3][0]) == r[3][0])
assert(pifrac*g*f == pifrac) #check that g*f is the identity on the function field of Y
"""
For convenience we work in the ring Pz with variables z0, ..., z8, which we identify with the rij
"""
Pz = PolynomialRing(QQ, ['z%s' %pasdf for pasdf in range(0, 9)], order='degrevlex')
Pz.inject_variables()
z = Pz.gens()[:9]
Pa = PolynomialRing(QQ, ['a%s' %pasdf for pasdf in range(0, 9)]+['t0', 't1','t2' ], order='degrevlex')
Pa.inject_variables()
a = Pa.gens()[:9]
QtoPz = Q.hom([z0,z1,z2,z3,z4,z5,z6,z7,z8], Pz)
QtoPzFrac = FractionField(Q).hom([z0,z1,z2,z3,z4,z5,z6,z7,z8], FractionField(Pz))
PztoQ = Pz.hom([r10, r20, r30, r11, r21, r31, r12, r22, r32 ], FractionField(Q))
"""
change all the rij's to zi's
"""
Yideal = QtoPz(Yideal)
eqlist = Yideal.gens()
"""
We now compute the subvariety of Y which is the inverse image under g of the points where two of the points 0,1,p0,p1,p2,q0,q1,q2, or where one of the points goes to infinity. We call an ideal cutting out a compontent representing the collision of two of these points a peripheral ideal.
"""
def get_peripheral_ideals():
gens = P.gens()
peripheral_ideals = set([])
for x in gens: #First find where a point collides with 0
poly = g(x).numerator()
factors = [x[0] for x in poly.factor()]
for factor in factors:
peripheral_ideals.add(QtoPz( Q.ideal(factor)))
for x in gens: #now where a point goes to infinity, i.e. one component of the denominator of g goes to 0
poly = g(x).denominator()
factors = [x[0] for x in poly.factor()]
for factor in factors:
peripheral_ideals.add(QtoPz(Q.ideal(factor)))
for x in gens: #now where a point collides with 1
poly = g(x-1).numerator()
factors = [x[0] for x in poly.factor()]
for factor in factors:
peripheral_ideals.add(QtoPz(Q.ideal(factor)))
for i in range(len(gens)): #now where two moving points collide
for j in range(i+1, len(gens)):
poly = g(gens[i] - gens[j] ).numerator()
factors = [x[0] for x in poly.factor()]
for factor in factors:
peripheral_ideals.add(QtoPz( Q.ideal(factor)))
return peripheral_ideals
"""
returns list of polynomials generating the peripheral ideals returned by get_peripheral_ideals
"""
def get_peripheral_polys():
ideals_tmp = get_peripheral_ideals()
## sage 6.6 bug workaround
ideals = []
for I in ideals_tmp :
if I not in ideals :
ideals.append(I)
##
polys = []
for I in ideals:
assert(len(I.gens()) == 1)
P = I.gens()[0]
polys.append(P)
return polys
try: #build peripheral polynomial list if it doesn't exist already
[periph_id, periph_poly] = load(file+'periph')
except IOError:
print 'building peripheral polynomials'
periph_id = get_peripheral_ideals()
periph_poly = get_peripheral_polys()
save([periph_id, periph_poly], file+'periph')
assert(len(periph_poly) == 35)
#check that the peripheral divisor has 35 components
"""
Build ResP and ResQ
ResP[i] is the residue of the 1-form at pi.
Ditto for ResQ
"""
try:
[ResP, ResQ] = load(file + 'Res')
except IOError:
print 'building residue equations'
ResP = [QtoPzFrac(g(p[j]*(p[j]-1)*(p[j] - y3)/((p[j] - q[j]) * prod([(p[j]-p[i]) * (p[j]-q[i]) for i in range(3) if i != j])))) for j in range(3)]
ResQ = [QtoPzFrac(g(q[j]*(q[j]-1)*(q[j] - y3)/((q[j] - p[j]) * prod([(q[j]-p[i]) * (q[j]-q[i]) for i in range(3) if i != j])))) for j in range(3)]
save([ResP, ResQ], file + 'Res')
def divide_by_peripheral(pol, periph):
"""
divide_by_peripheral: takes polynomials pol and periph returns pol after dividing by periph
until it divides no more
"""
ring = pol.parent()
for J in periph:
while(J.divides(pol)):
pol = ring(pol/J)
return pol
def build_OppRes():
"""
build the two polynomials expressing the opposite reside condition
"""
opp = [(ResP[j]+ResQ[j]).numerator() for j in [0,1]]
opp2 = [divide_by_peripheral(pol, periph_poly) for pol in opp]
return opp2
"""
Build the opposite residue equations.
These equations express the condition that a point in Y' represents a form which satisfies the opposite-residue condition
"""
try: #build opposite residue equations list if it doesn't exist already
OppRes = load(file+'OppRes')
except IOError:
print 'building opposite residue equations'
OppRes = build_OppRes()
save(OppRes, file+'OppRes')
"""
Build the list of subspaces constructed by the first part of the torus containment algorithm.
These are the potential subspaces which could yield a torus-translate contained inY t
"""
try: #build subspaces list if it doesn't exist already
subspaces = load(file+'subspaces')
except IOError:
print 'building subspaces'
subspaces = build_all_subspaces(matrix_list, eqlist)
save(subspaces, file+'subspaces')
print 'rank 3 subspaces:', len(subspaces[0])
print 'rank 2 subspaces:', len(subspaces[1])
print 'rank 1 subspaces:', len(subspaces[2])
def has_lonely(eqlist, N):
"""
returns True if the partition of the exponent vectors of any polynomial in eqlist has a part
consisting of a single vector. This rules out any nonperipheral torus-translate associated to
N.
"""
from collections import Counter
for eq in eqlist:
projection_vectors = [vector(x) for x in eq.exponents()]
projections = [orthogonal_projection(v, N) for v in projection_vectors]
for v in projections:
v.set_immutable()
count = Counter(projections).values()
if 1 in count:
return True
subspaces_to_handle = flatten(subspaces)
print 'total subspaces:', len(subspaces_to_handle)
try:
handled_by_lonely = load(file + 'handled_by_lonely')
except IOError:
handled_by_lonely = []
for N in subspaces_to_handle:
if has_lonely(eqlist, N):
handled_by_lonely.append(N)
save(handled_by_lonely, file + 'handled_by_lonely')
def diff(a, b):
"""
elements of list a not in list b
"""
b = set(b)
return [aa for aa in a if aa not in b]
subspaces_to_handle = diff(subspaces_to_handle, handled_by_lonely)
print len(handled_by_lonely), ' subspaces handled by lonely vector condition.'
print len(subspaces_to_handle), ' subspaces remain'
print ''
def torus_hom(N):
"""
ring homomorphism induced by the parameterized torus-translate (a, t) |--> a t^N
"""
if len(N.rows()) == 1:
e = N.rows()[0]
f = (ZZ^9).zero_vector()
g = (ZZ^9).zero_vector()
elif len(N.rows()) == 2:
e = N.rows()[0]
f = N.rows()[1]
g = (ZZ^9).zero_vector()
elif len(N.rows()) == 3:
e = N.rows()[0]
f = N.rows()[1]
g = N.rows()[2]
else:
raise Exception
image = [a[i]*t0^e[i] * t1^f[i]*t2^g[i] for i in range(9)]
return FractionField(Pz).hom(image, FractionField(Pa))
def coefficient_ideal(equations, N):
"""
equations are a list of equations in the zi
N the matrix defining a torus in the zi
returns ideal in Q[ai] yielding torus translates contained in the variety cut out by equations
"""
assert len(N.rows()) in range(1,4) #only implemented for these cases
TH = torus_hom(N)
"""
find the list of polynomials, new_eqs, which cut out TH^{-1}(Y) where Y is cut out by equations.
i.e. we find polynomials cutting out those (a, t) such that a t^N \in Y
"""
new_eqs = [TH(eq).numerator() for eq in equations]
new_eqs2 = []
"""
We wish to find a such that a t^N \in Y for all t. These are exactly the a such that p(a,t)=0
for all t and for all p in new_eqs. We consider each polynomial in new_eqs as an element of
Q[a][t]. For each monomial in t, the condition is that the coefficient, a poly in Q[a],
vanishes. So we collect below these coefficients.
"""
for new_eq in new_eqs:
deg = new_eq.degree()
for d0 in range(deg+1):
for d1 in range(deg + 1):
for d2 in range(deg + 1):
tempeq = new_eq.coefficient({t0:d0, t1:d1, t2:d2}) #coefficient of t0^d0 t1^d1 t2^d2
if tempeq != 0:
new_eqs2.append(tempeq)
return Pa.ideal(new_eqs2)
def peripheral_quotients(eqlist, N, verbose = False):
"""
eqlist are a list of equations in the zi
N the matrix defining a torus in the zi
returns ideal in Q[ai] yielding torus translates contained in the variety cut out by equations,
after removing irreducible components which yield peripheral tori (in the boundary of M_0,10)
"""
if verbose:
print ''
print N
"""
I is the ideal cutting out the variety of coefficients a such that the torus-translate a*t^N is
contained in the variety cut out by eqlist.
"""
I = coefficient_ideal(eqlist, N)
TH = torus_hom(N)
"""
list of polynomials cutting out inverse image of peripheral divisors under (a,t) |-->at^N, as
well as coordinate axes.
"""
periph_ideal_a = [Pa.ideal(TH(poly).numerator()) for poly in periph_poly] + [Pa.ideal(a[i]) for i in range(9)]
"""
remove irreducible components which parameterize peripheral tori
"""
for J in periph_ideal_a:
I = I.saturation(J)[0]
return I
def peripheral_quotients_in_Pa(I, N):
"""
similar to peripheral_quotients. Here I is an ideal in Q(ai) (rather than pulling back an ideal in Q(zi) in peripheral_quotients).
Returns the saturation ideal I:J where J is obtained by pulling back the peripheral ideal in Q(zi).
In words, we are finding the open subset of V(I) consisting of a such that the torus-translate
a*t^N is nonperipheral.
"""
TH = torus_hom(N)
periph_ideal_a = [Pa.ideal(TH(poly).numerator()) for poly in periph_poly] + [Pa.ideal(a[i]) for i in range(9)]
for J in periph_ideal_a:
I = I.saturation(J)[0]
return I
def Opp_Res_ideal(N):
"""
Ideal expressing condition on coefficient (a) that the torus at^N is contained in the locus of
forms satisfying the opposite residue condition
"""
return coefficient_ideal(OppRes, N)
def Const_Res_ideal(N):
assert(len(N.rows()) == 1)
TH = torus_hom(N)
rat = TH(ResP[0]/ResP[1])
const_rat = rat.derivative(t0).numerator()
deg = const_rat.degree()
crI1 = Pa.ideal([const_rat.coefficient({t0:d}) for d in range(deg)])
rat = TH(ResP[0]/ResP[2])
const_rat = rat.derivative(t0).numerator()
deg = const_rat.degree()
crI2 = Pa.ideal([const_rat.coefficient({t0:d}) for d in range(deg)])
return crI1 + crI2
def build_basic_ideals(subspaces_to_handle, eqlist):
basic_ideals = {}
for N in subspaces_to_handle:
basic_ideals[N] = peripheral_quotients(eqlist, N)
return basic_ideals
try: #build dictionary of ideals using only equations cutting out M_0,10 if it doesn't exist already
basic_ideals = load(file+'basic_ideals')
except IOError:
print 'building basic ideals'
basic_ideals = build_basic_ideals(subspaces_to_handle, eqlist)
save(basic_ideals, file + 'basic_ideals')
try:
handled_by_basic_ideals = load(file + 'handled_by_basic_ideals')
except IOError:
handled_by_basic_ideals = []
for N in subspaces_to_handle:
if basic_ideals[N] == Pa.ideal(1):
handled_by_basic_ideals.append(N)
save(handled_by_basic_ideals, file + 'handled_by_basic_ideals')
print len(handled_by_basic_ideals), 'subspaces handled by saturation ideals'
subspaces_to_handle = diff(subspaces_to_handle, handled_by_basic_ideals)
print len(subspaces_to_handle), 'subspaces remain'
print ''
try:
handled_by_opp_res = load(file + 'handled_by_opp_res')
opp_res_ideals = load(file + 'opp_res_ideals')
except IOError:
handled_by_opp_res = []
opp_res_ideals = {}
for N in subspaces_to_handle:
if N != matrix([[0, 1, 1, 0, 1, 1, 0, 1, 1]]): #this one is handled separately below
I = basic_ideals[N]
J = Opp_Res_ideal(N)
L = I + J
M = peripheral_quotients_in_Pa(L, N)
if M == Pa.ideal(1):
handled_by_opp_res.append(N)
else:
opp_res_ideals[N] = M
save(handled_by_opp_res, file + 'handled_by_opp_res')
save(opp_res_ideals, file + 'opp_res_ideals')
subspaces_to_handle = diff(subspaces_to_handle, handled_by_opp_res)
print len(handled_by_opp_res), 'handled by opposite residue coondition'
print len(subspaces_to_handle), 'subspaces remain'
"""
Now consider (0,1,1,0,1,1,0,1,1).
This case is alse handled using just the M_0,10 equations together with the opposite residue equations.
For some reason the computation only completes if we compute the primary decomposition of the basic_ideal first, then apply the previous strategy to each associated prime, so this case is handled separately
"""
N = matrix([[0,1,1,0,1,1,0,1,1]])
N.set_immutable()
assert(N in subspaces_to_handle)
primary = basic_ideals[N].primary_decomposition_complete()
primes = [x[1] for x in primary] #list of associated primes
assert(len(primes) == 6)
J = Opp_Res_ideal(N)
for I in primes:
L = I+J
assert(peripheral_quotients_in_Pa(L, N) == Pa.ideal(1))
subspaces_to_handle.remove(N)
print "handled special case:", N
assert(subspaces_to_handle == []) #We're done!!!
print("Wall time: "+str(walltime()-start_time)+"s")
print ':D' #yay