MO442252: enumeration of Deligne-Mostow lattices

orbifold.py

#!/usr/bin/pypy3

"""
Ethan Dlugie asked [1] "How do we know there are no more Deligne-Mostow/Thurston lattices?"
than the 94 enumerated by Thurston [3]. This program combines case analysis with exhaustive search
to show that the list is in fact complete.

There is a subtle difference between Thurston's claim (Theorem 0.2) and Mostow's ΣINT condition [2]:
Mostow only allows one set of equal angles to give half-integer reciprocals (condition (c) below),
whereas Thurston allows multiple sets. This code uses Thurston's condition, as more straightforward,
but normalised following Mostow, again for simplicity.

Then we want to enumerate partitions of 2 into five or more real numbers mu_i in (0, 1) for which
every pair of indices i != j satisfies one of three conditions:
 (a) mu_i + mu_j >= 1; or
 (b) 1 / (1 - mu_i - mu_j) is an integer; or
 (c) mu_i = mu_j and 2 / (1 - mu_i - mu_j) is an integer.

Dlugie noted that a possible strategy would be to generate partitions into five parts and then to
refine them to get partitions into six parts, etc. This is the strategy followed herein.

On notation: define r_{i,j} as
	r_{i,j} = 1 / (1 - mu_i - mu_j)
Inverting this we get
	mu_i + mu_j = 1 - 1 / r_{i,j}
which will be used heavily. Note that r_{i,j} will only be used in contexts where it's a positive half-integer.

Since we're dealing with partitions we can assume without loss of generality that they're in descending order.
This assumption is useful not only for avoiding the need to deduplicate permuted solutions but also for bounding
undetermined parts in the case analyses.


[1] https://mathoverflow.net/q/442252
[2] Mostow, G. D., Generalized Picard lattices arising from half-integral conditions,
    Publ. Math., Inst. Hautes Étud. Sci. 63, 91-106 (1986). ZBL0615.22009.
[3] Thurston, William P., Shapes of polyhedra and triangulations of the sphere, Rivin, Igor (ed.) et al.,
    The Epstein Birthday Schrift dedicated to David Epstein on the occasion of his 60th birthday.
	Geom. Topol. Monogr. 1, 511-549 (1998). ZBL0931.57010.
"""

from fractions import Fraction
from itertools import product
from math import gcd

def validate(soln):
	"""
	Verifies that a candidate solution satisfies the conditions: it should have the correct sum,
	the parts should be in the correct interval, and the pairs of parts should meet one of the three
	conditions (a), (b), (c).
	"""
	if sum(soln) != 2 or min(soln) <= 0 or max(soln) >= 1:
		return False

	while soln:
		mu_i = soln[0]
		soln = soln[1:]
		for mu_j in soln:
			if mu_i + mu_j < 1:
				r = rbound(mu_i, mu_j)
				if not (r == int(r) or mu_i == mu_j and 2*r == int(2*r)):
					return False

	return True


def rbound(a, b):
	""" Given upper bounds on mu_i and mu_j, returns an upper bound on r_{i,j} """
	assert a + b < 1
	return Fraction(1) / (1 - a - b)


def rrange(rbound, rbound2 = None):
	"""
	Enumerates possible values of r_{a,b} in a range.
	Note that unlike the Python built-in range, this takes its upper bound as inclusive.
	Can be called with just an upper bound or with a lower bound and an upper bound.
	"""
	# Decode the arguments. Since we're dealing with half-integers, we want to double and round appropriately.
	if rbound2 is None:
		# Lower bound is 3/2, the smallest sensible value of r
		lo = 3
		# Upper bound is rbound rounded down to the nearest half-integer.
		hi = int(2*rbound)
	else:
		# Lower bound is rbound rounded up to the nearest half-integer.
		lo = 2*rbound
		lo = int(lo) + (0 if lo == int(lo) else 1)
		assert lo >= 3

		# Upper bound is rbound2 rounded down to the nearest half-integer.
		hi = int(2*rbound2)

	for n in range(lo, hi + 1):
		yield Fraction(n, 2)


def mu_sum(r):
	""" Takes r_{a,b}, which must be a half-integer greater than 1, and returns mu_a + mu_b """
	assert r > 1
	assert 2*r == int(2*r), str(r)
	return 1 - Fraction(1, r)


def quins345(mu_3, mu_4, mu_5):
	""" Given mu_3, mu_4, mu_5 generates all completions to a solution in five parts. """
	assert mu_3 >= mu_4 >= mu_5
	# Since mu_1 >= mu_2, mu_2 is at most half of the unallocated total weight.
	mu2_bound = (2 - mu_3 - mu_4 - mu_5) / 2
	# This assertion must hold because otherwise mu2_bound + mu_4 is also at least 1, so 2 mu2_bound + mu_3 + mu_4 + mu_5 > 2.
	assert mu2_bound + mu_5 < 1
	for r25 in rrange(rbound(mu2_bound, mu_5)):
		mu_2 = mu_sum(r25) - mu_5
		mu_1 = 2 - mu_2 - mu_3 - mu_4 - mu_5
		candidate = (mu_1, mu_2, mu_3, mu_4, mu_5)
		if mu_1 >= mu_2 >= mu_3 and validate(candidate):
			yield candidate


def mus_from_rs(r_ij, r_ik, r_jk):
	""" Given the r values for all pairs of a triple, returns the corresponding mu values """
	assert r_ij >= r_ik >= r_jk
	sm_ij = mu_sum(r_ij)
	sm_ik = mu_sum(r_ik)
	sm_jk = mu_sum(r_jk)
	sm_ijk = (sm_ij + sm_ik + sm_jk) / 2
	return sm_ijk - sm_jk, sm_ijk - sm_ik, sm_ijk - sm_ij


def quins():
	""" Enumerates all solutions in five parts """
	# We assume mu_1 >= mu_2 >= ... >= mu_5. Then mu_4 + mu_5 <= 4/5 so r_{4,5} <= 5.
	for r45 in rrange(5):
		sm45 = mu_sum(r45)
		# mu_3 is at most a third of the remaining total
		mu3_bound = (2 - sm45) / 3
		if r45 != int(r45):
			# Half-integer, so mu_4 = mu_5
			mu_4 = sm45 / 2
			mu_5 = mu_4
			for r35 in rrange(r45, rbound(mu3_bound, mu_5)):
				mu_3 = mu_sum(r35) - mu_5
				assert mu_3 >= mu_4
				yield from quins345(mu_3, mu_4, mu_5)
		elif r45 > 2:
			# In this case the bound on mu_3 is low enough to make things easy.
			assert mu3_bound < Fraction(1, 2)
			# mu_4 >= mu_5 so mu_5 is at most half the sum
			for r35 in rrange(r45, rbound(mu3_bound, sm45 / 2)):
				# mu_3 >= mu_4 so mu3_bound bounds both
				for r34 in rrange(r35, rbound(mu3_bound, mu3_bound)):
					yield from quins345(*mus_from_rs(r34, r35, r45))
		else:
			# One case to cover.
			assert r45 == 2
			assert sm45 == Fraction(1, 2)
			assert mu3_bound == Fraction(1, 2)

			# mu_5 is at most half of mu_4 + mu_5, and we bound r_{3,5} <= 4.
			r35_bound = rbound(mu3_bound, sm45 / 2)
			assert r35_bound == 4

			# Case r_{3,5} is half-integer but not integer: then mu_3 = mu_5.
			# Since they're in descending order mu_3 = mu_4 = mu_5.
			# But then r_{3,5} = r_{4,5} = 2 giving a contradiction.

			# Case r_{3,5} = 4: then both mu_3 and mu_5 are at their bounds, so mu_4 = mu_5 = 1/4 and mu_3 = 1/2.
			# That leaves total weight 1 to split between mu_1 and mu_2 such that mu_1 >= mu_2 >= mu_3, forcing all values.
			assert mu_sum(4) == Fraction(3, 4)
			candidate = (Fraction(1, 2), Fraction(1, 2), Fraction(1, 2), Fraction(1, 4), Fraction(1, 4))
			assert validate(candidate)
			yield candidate

			# Case r_{3,5} = 3: the bound on mu_3 gives an adequate bound on mu_4
			r35 = 3
			mu4_bound = sm45 - (mu_sum(r35) - mu3_bound)
			for r34 in rrange(r35, rbound(mu3_bound, mu4_bound)):
				yield from quins345(*mus_from_rs(r34, r35, r45))

			# Case r_{3,5} = 2: this is the tricky subcase.
			r35 = 2
			# Since r_{3,5} = r_{4,5} we have mu_3 = mu_4. But mu_4 >= sm45/2 so we can bound mu_2 <= 5/8.
			mu2_bound = (2 - Fraction(3, 2) * sm45) / 2
			assert mu2_bound == Fraction(5, 8)
			for r25 in rrange(r35, rbound(mu2_bound, sm45 / 2)):
				# r_{2,5} must be integral for the same reason that r_{3,5} must be integral.
				if r25 != int(r25):
					continue

				# We require mu_3 = mu_4 < 1/2 because otherwise we get mu_5 = 0.
				# So mu_3 + mu_4 < 1 and 2 r_{3,4} is an integer, which we shall denote k.
				# From the definition of r_{a,b} and the equality of mu_3 and mu_4 we get mu_4 = 1/2 - 1/k.
				# But we also have mu_4 + mu_5 = 1/2, so mu_5 = 1/k. Note that k >= 4 since mu_4 >= mu_5.

				if r25 == 2:
					# mu_2 = mu_4, and by the known total we have mu_1 = 1/2 + 2/k.
					# If mu_1 + mu_5 = 1/2 + 3/k < 1 then r_{1,5} = 2 + 12/(k-6) is an integer.
					# So the candidates are 4, 5, 6 (for which mu_1 + mu_5 >= 1) and those
					# 7 <= k <= 18 for which (k-6) | 12.
					for k in range(4, 19):
						candidate = (Fraction(k+4,2*k), Fraction(k-2,2*k), Fraction(k-2,2*k), Fraction(k-2,2*k), Fraction(2,2*k))
						if validate(candidate):
							yield candidate
				else:
					#   mu_4 = mu_3 = 1/2 - 1/k
					#   mu_2 + mu_5 = 1 - 1 / r_{2,5}
					#   mu_1 = 2 - mu_2 - mu_3 - mu_4 - mu_5 >= mu_2
					# So
					#   2/k + 1 / r_{2,5} >= mu_2
					# Add mu_5 = 1/k to each side to get
					#   3/k + 1 / r_{2,5} >= 1 - 1 / r_{2,5}
					# and rearrange to
					#   3/k >= 1 - 2 / r_{2,5}
					# so that
					#   k <= 3r_{2,5} / (r_{2,5} - 2)
					for k in range(4, int(3*r25 /(r25 - 2)) + 1):
						mu_5 = Fraction(1, k)
						mu_4 = sm45 - mu_5
						mu_3 = mu_4
						mu_2 = mu_sum(r25) - mu_5
						mu_1 = 2 - mu_2 - mu_3 - mu_4 - mu_5
						candidate = (mu_1, mu_2, mu_3, mu_4, mu_5)
						if validate(candidate):
							yield candidate


def orbifolds():
	""" Generates all orbifold partitions of 2 into five or more parts """
	# Helper: putting solutions into canonical order avoids duplication.
	def canonical(soln):
		return tuple(reversed(sorted(soln)))

	# Start by generating the quintuplet solutions
	layer = set(quins())
	while layer:
		# Emit the solutions for the current level of refinement
		yield from layer

		# Acculumator for refined solutions with one more part. Use a set for deduplication.
		# Note that we accumulate without testing for validity. We'll filter later.
		next_layer = set()
		for soln in layer:
			# Pick a part to split
			for i, mu_i in enumerate(soln):
				# The parts other than the one we split
				others = soln[:i] + soln[i+1:]

				# If we split mu_i, we get two parts which sum to mu_i
				rsplit = 1 / (1 - mu_i)
				# Bound the smaller part of the split
				h = mu_i / 2
				if rsplit == int(rsplit):
					# Find the largest unsplit part which gives a sum less than unity,
					# because that will give us the tightest possible bound on r.
					mu_k = max(m for m in others if m + h < 1)
					for r in rrange(rbound(mu_k, h)):
						split_lo = mu_sum(r) - mu_k
						next_layer.add(canonical(others + (split_lo, mu_i - split_lo)))
				elif 2*rsplit == int(2*rsplit):
					# Has to split into equal parts
					next_layer.add(canonical(others + (h, h)))

		# Now filter for validity.
		layer = [candidate for candidate in next_layer if validate(candidate)]


def format_for_printing(soln):
	""" Converts the fractions into the common denominator and a list of numerators, as in Thurston's appendix """
	denom = 1
	for mu in soln:
		denom *= mu.denominator // gcd(denom, mu.denominator)

	return denom, tuple(int(mu * denom) for mu in soln)


if __name__ == "__main__":
	for i, (denom, numerators) in enumerate(sorted(format_for_printing(soln) for soln in orbifolds())):
		print(i+1, "\t", denom, "\t", numerators)