from copy import deepcopy from math import sqrt import sys _verbose = True ############################################################################### # Section 1 # # The tetrahedron class class Tetrahedron: ''' Represents a tetrahedron of a triangulation where each face gluing is by the identity, i.e., vertex i of face j with i != j will be glued to vertex i of face j of another tetrahedron. Thus the gluing is determined by only one pointer to another tetrahedron per face. Besides the four pointers to the glued tetrahedra, each tetrahedron carries an orientation which can be positive or negative and which is set when a tetrahedron is instantiated with Tetrahedron(orientation = +1), respectively, Tetrahedron(orientation = -1) and can be read by tet.GetOrientation(). One basic operations to do for instances of this tetrahedra class is GlueFaces(A, B, f) which just glues tetrahedron A to tetrahedron B along face f. It is required that A and B have opposite orientations. The other basic operation is a feature not supported by SnapPy or Regina triangulation: Identify(A, B) will identify two tetrahedra A and B and the identification is pushed through the face pairings. E.g., if A is glued to A' along face 0 and B to B' along face 0, then A' and B' also will be identified and similarly for other tets glued to A' and B' and so on. This is useful for constructing, e.g., a quotient space by a group action: one can fix a tetrahedron of the triangulation, compute its images under the group action and then identify the tetrahedron with all of its images. Note that the tetrahedra A and B above will be the same tetrahedron in the quotient though there will still be different as python objects. Thus, one has to use the API AreIdentified(tet1, tet2) to determine whether two tetrahedra are the same or not after an identification. So far, all API was concerned about setting/constructing triangulations. The API to "read" triangulations is tet.GetFaceGluing(f) which returns the tetrahedron glued to tet along face f. Based on this API, there are further convenience functions such as tet.Rotate(t, faceA, faceB), tet.Order(faceA, faceB) and tet.Orders(). Assume a client had two references to tetrahedra A and B and they get identified, either directly by an Identify(A, B) call or because some of other tetrahedra got identified that cause that A and B need to be identified. The data structures around the tetrahedra will always maintain a consistent triangulation. For example, AreIdentified(A.GetFaceGluing(0).GetFaceGluing(0)) will always be True unless face 0 of A was never glued. Furthermore, as A and B are identified, all public API calls on A and B will return the same answers. The user should think of them as practically the same, just different representations as python references. This is, however, not true for the "is" and "==" operator, so the user should always use AreIdentified instead of A == A.GetFaceGluing(0).GetFaceGluing(0) . If the user is holding on to a bunch of references and several might actually redundantly point to the same identified tetrahedron, he or she might wish to consolidate the set/list into a set that only contains a unique reference per identified tetrahedron. This is done by the FilterCanonical function. I.e., FilterCanonical([A,B]) for the tetrahedra above will returns a set of exactly one python object representing the identified tetrahedra that both A and B were representing. The result might not be the python reference stored in A or B if A and B were at some identified with yet another tetrahedron. If the user created a huge set S of references to tetrahedra that have been identified, he or she is encouraged to run S = FilterCanonical(S). This will not change the set of identified tetrahedra, but it will facilitate garbage collecting redundant references to those tetrahedra. >>> F1 = Tetrahedron(orientation = +1) >>> F2 = Tetrahedron(orientation = -1) >>> GlueFaces(F1, F2, 0) True >>> F3 = Tetrahedron(orientation = +1) >>> Identify(F1, F3) True >>> AreIdentified(F1, F3) True >>> AreIdentified(F3.GetFaceGluing(0), F2) True >>> CheckGluingConsistency(GetConnectedComponent(F3)) ''' def __init__(self, orientation): self._faceGluings = 4 * [ None ] self._orientation = orientation self._canonicalRepresentative = self def GetOrientation(self): return self._orientation def GetFaceGluing(self, face): f = self.GetCanonicalRepresentative()._faceGluings[face] if f: self._faceGluings[face] = f.GetCanonicalRepresentative() return self._faceGluings[face] def GetCanonicalRepresentative(self): if self is self._canonicalRepresentative: return self self._canonicalRepresentative = ( self._canonicalRepresentative.GetCanonicalRepresentative()) return self._canonicalRepresentative def Rotate(self, t, faceA, faceB): if t > 0: return (self.Rotate(t-1, faceA, faceB) .GetFaceGluing(faceA).GetFaceGluing(faceB)) if t < 0: return self.Rotate(-t, faceB, faceA) return self def Order(self, faceA, faceB): try: r = 1 tet = self.Rotate(1, faceA, faceB) while not self.GetCanonicalRepresentative() is tet: r += 1 tet = tet.Rotate(1, faceA, faceB) return r except: return -1 def Orders(self): return (self.Order(0,1), self.Order(1,2), self.Order(2,3), self.Order(0,2), self.Order(0,3), self.Order(1,3)) def EnforceEdgeCycle(self, degree, faceA, faceB): tet = self for i in range(2 * degree + 2): nextTet = tet.GetFaceGluing(faceA if i % 2 == 0 else faceB) if not nextTet: break else: tet = nextTet return ( tet._EnforceEdgeCycle(degree, faceA, faceB) or tet._EnforceEdgeCycle(degree, faceB, faceA)) def _EnforceEdgeCycle(self, degree, faceA, faceB): tet = self for i in range(2 * degree - 1): tet = tet.GetFaceGluing(faceA if i % 2 == 0 else faceB) if not tet: return False next = tet.GetFaceGluing(faceB) prev = self.GetFaceGluing(faceB) if next: return Identify(self, next) elif prev: return Identify(prev, tet) else: return GlueFaces(tet, self, faceB) ############################################################################### # Section 2 # # Basic functions for tetrahedra def GlueFaces(tet1, tet2, face): """ Glues to tetrahedra tet1 and tet2 along face with index face. """ # Get canonical representatives for tetrahedra canonTet1 = tet1.GetCanonicalRepresentative() canonTet2 = tet2.GetCanonicalRepresentative() # Make sure they were not glued before assert not canonTet1._faceGluings[face] assert not canonTet2._faceGluings[face] # Make sure the orientations are compatible assert canonTet1.GetOrientation() == -canonTet2.GetOrientation() # Glue canonTet1._faceGluings[face] = canonTet2 canonTet2._faceGluings[face] = canonTet1 return True def Identify(tet1, tet2): """ Identifies tetrahedron tet1 and tet2. Returns true if and only if the two tetrahedra had not been identified before. """ # Get canonical representatives for tetrahedra canonTet1 = tet1.GetCanonicalRepresentative() canonTet2 = tet2.GetCanonicalRepresentative() # Make sure the orientations are compatible assert canonTet1.GetOrientation() == canonTet2.GetOrientation() if canonTet1 is canonTet2: # Do nothing when tets were already identified return False # A working set of pairs (tet1, tet2) that need to be identified # Start with the given pair, then push more pairs onto it when # tetrahedra that are identified were glued to other tetrahedra pairs = [(canonTet1, canonTet2)] while pairs: # Get the next pair of tetrahedra to be idenfied tet1, tet2 = pairs.pop() canonTet1 = tet1.GetCanonicalRepresentative() canonTet2 = tet2.GetCanonicalRepresentative() # Move on to next pair if alreay identified if canonTet1 is canonTet2: continue # Identify: this means that the second tetrahedron sets # its pointer to the canonical representative of the # equivalence class of identified tetrahedra. It sets this # pointer to the first tetrahedron canonTet2._canonicalRepresentative = canonTet1 # Push the identification through the faces for face in range(4): # Get old face gluings of tets face1 = canonTet1._faceGluings[face] if face1: face1 = face1.GetCanonicalRepresentative() face2 = canonTet2._faceGluings[face] if face2: face2 = face2.GetCanonicalRepresentative() if not face1: # If the first tetrahedron was not glued but the # second tetrahedron was glued to face2, we # glue the first tetrahedron to face2. # The first tetrahedron becomes the canonical representative # so the value of _faceGluings of the second tetrahedron # will never be accessed again. That's why we do not care # setting it. canonTet1._faceGluings[face] = face2 elif face2: # If both tets were glued to tetrahedra, we need to identify # the glued tetrahedra as well. Push them onto the working set. pairs.append((face1, face2)) return True def AreIdentified(tet1, tet2): """ Returns whether tet1 and tet2 were identified """ return ( tet1.GetCanonicalRepresentative() is tet2.GetCanonicalRepresentative()) ############################################################################### # Section 3 # # Utility functions for tetrahedra def FilterCanonical(tets): """ Given a set of tetrahedra, return a set containing exactly one representative for each equivalence class of identified tetrahedra. """ return set([tet.GetCanonicalRepresentative() for tet in tets]) def FirstPositivelyOrientedTet(tets): """ Get the first positively oriented tetrahedra from the iterable tets. For a list, this is deterministically the first one. For a set, it is any tetrahedron. """ for tet in tets: if tet.GetOrientation() > 0: return tet def FirstNegativelyOrientedTet(tets): """ Similar to FirstPositivelyOrientedTet, but opposite orientation. """ for tet in tets: if tet.GetOrientation() < 0: return tet def GetConnectedComponent(tet, faces = [0, 1, 2, 3]): """ Get the connected component of a triangulation containing tet. This recursively follows all faces of the given tetrahedron. If faces are given, it only follows the given faces. """ tets = set() workingSet = [ tet.GetCanonicalRepresentative() ] while workingSet: tet = workingSet.pop() if not tet in tets: tets.add(tet) for face in faces: otherTet = tet.GetFaceGluing(face) if otherTet: workingSet.append(otherTet) return tets def CheckGluingConsistency(tets): """ Checks that the triangulation formed by tets is consistent. If only the public API is used, triangulations should always be in a consistent state. """ for tet in tets: for face in range(4): g1 = tet.GetFaceGluing(face) if g1: g2 = g1.GetFaceGluing(face) assert not g1 is g2 assert g2 is tet.GetCanonicalRepresentative() def TetrahedraToReginaTriangulation(tets): """ Given a set of a tetrahedra forming a triangulation, create a regina triangulation from it. open('SnapPeaTriangulation.trig' % (d, i),'w').write( TetrahedronToReginaTriangulation(tets).snapPea()) """ import regina tets = FilterCanonical(tets) reginaTriangulation = regina.NTriangulation() # Create a regina tetrahedron for each given tetrahedron tetsToReginaTets = dict([ (tet, reginaTriangulation.newTetrahedron()) for tet in tets]) for tet in tets: for face in range(4): # Glue the regina tetrahedra using trivial face permutation otherTet = tet.GetFaceGluing(face) if otherTet: assert tetsToReginaTets.has_key(otherTet), ( "Supplied set of tetrahedra is not closed under " "face gluings.") reginaTet = tetsToReginaTets[tet] otherReginaTet = tetsToReginaTets[otherTet] reginaTet.joinTo(face, otherReginaTet, regina.NPerm4()) return reginaTriangulation ############################################################################### # Section 4 # # Arithmetic classes and helpers class AlgebraicInteger4(object): """ Represents an element in the ring of integers of the imaginary quadratic number field of discriminant -4. This number is hashable and comparable so that it can be used in sets and to find a canonical element among the set with "min". Notice that the comparison is deterministic but not natural as there is no natural ordering on this ring. """ def __init__(self, a, b): """ The number a + bi """ self.a = a self.b = b def Conjugate(self): return AlgebraicInteger4(self.a, -self.b) def NormSquared(self): return (self * self.Conjugate()).a def __add__(self, other): return AlgebraicInteger4(self.a + other.a, self.b + other.b) def __neg__(self): return AlgebraicInteger4(-self.a, -self.b) def __sub__(self, other): return self + (-other) def __mul__(self, other): return AlgebraicInteger4( self.a * other.a - self.b * other.b, self.a * other.b + self.b * other.a) def __mod__(self, other): """ a % b returns a canonical element in the congruence class mod b. """ conjugate = other.Conjugate() normSquared = (other * conjugate).a prod = self * conjugate m = AlgebraicInteger4(prod.a / normSquared, prod.b /normSquared) return self + (-m * other) def __eq__(self, other): if isinstance(other, int): return self.a == other and self.b == 0 return self.a == other.a and self.b == other.b def __hash__(self): return hash((self.a, self.b)) def __cmp__(self, other): return cmp((self.a, self.b), (other.a, other.b)) @staticmethod def Units(): return [ AlgebraicInteger4(1,0), AlgebraicInteger4(0,1), AlgebraicInteger4(0,-1), AlgebraicInteger4(-1,0) ] def __str__(self): return '%d + %d * I' % (self.a, self.b) class AlgebraicInteger3(object): """ Represents an element in the ring of integers of the imaginary quadratic number field of discriminant -3. This number is hashable and comparable so that it can be used in sets and to find a canonical element among the set with "min". Notice that the comparison is deterministic but not natural as there is no natural ordering on this ring. """ def __init__(self, a, b): self.a = a self.b = b def Conjugate(self): return AlgebraicInteger3(self.a + self.b, -self.b) def NormSquared(self): return (self * self.Conjugate()).a def __add__(self, other): return AlgebraicInteger3(self.a + other.a, self.b + other.b) def __neg__(self): return AlgebraicInteger3(-self.a, -self.b) def __sub__(self, other): return self + (-other) def __mul__(self, other): return AlgebraicInteger3( self.a * other.a - self.b * other.b, self.a * other.b + self.b * other.a + self.b * other.b) def __mod__(self, other): """ a % b returns a canonical element in the congruence class mod b. """ conjugate = other.Conjugate() normSquared = (other * conjugate).a prod = self * conjugate m = AlgebraicInteger3(prod.a / normSquared, prod.b /normSquared) return self + (-m * other) def __eq__(self, other): if isinstance(other, int): return self.a == other and self.b == 0 return self.a == other.a and self.b == other.b def __hash__(self): return hash((self.a, self.b)) def __cmp__(self, other): return cmp((self.a, self.b), (other.a, other.b)) @staticmethod def Units(): return [ AlgebraicInteger3(1,0), AlgebraicInteger3(0,1), AlgebraicInteger3(-1,1), AlgebraicInteger3(-1,0), AlgebraicInteger3(-1,0), AlgebraicInteger3(1,-1) ] def __str__(self): return '%d + %d * u' % (self.a, self.b) def UnitsInRingModZ(z): """ Get a representative for each unit in the ring O_d/z """ # Units in the ring of integers O_d units = z.Units() # Get all values in the ring O_d/z maxEntry = int(sqrt(z.NormSquared()) + 2) allValues = set([ type(z)(a, b) % z for a in range(-maxEntry, maxEntry) for b in range(0, 2 * maxEntry) ]) def MyMin(l): return min(l, key = lambda x: x.NormSquared()) def ShortestLinearCombination(a, b): return MyMin([ a - unit * b for unit in units ]) def Gcd(a, b): while True: if a == 0: return b if b == 0: return a a, b = MyMin([a, b]), ShortestLinearCombination(a, b) def IsUnit(a): return Gcd(z, a) in units return set([ a for a in allValues if IsUnit(a)]) class Matrix2x2: """ Represents a 2x2 matrix. Similar to AlgebraicIntegers3/4, objects are hashable and comparable. """ def __init__(self, z11, z12, z21, z22): self.z11 = z11 self.z12 = z12 self.z21 = z21 self.z22 = z22 def __mul__(self, other): """ Matrix multiplication or multiplication of matrix with a scalar. """ if isinstance(other, Matrix2x2): return Matrix2x2( self.z11 * other.z11 + self.z12 * other.z21, self.z11 * other.z12 + self.z12 * other.z22, self.z21 * other.z11 + self.z22 * other.z21, self.z21 * other.z12 + self.z22 * other.z22) else: return Matrix2x2( self.z11 * other, self.z12 * other, self.z21 * other, self.z22 * other) def __mod__(self, other): return Matrix2x2( self.z11 % other, self.z12 % other, self.z21 % other, self.z22 % other) def Det(self): return self.z11 * self.z22 - self.z12 * self.z21 def __str__(self): return '[ %s, %s; %s, %s]' % (self.z11, self.z12, self.z21, self.z22) def __hash__(self): return hash((self.z11, self.z12, self.z21, self.z22)) def __cmp__(self, other): return cmp((self.z11, self.z12, self.z21, self.z22), (other.z11, other.z12, other.z21, other.z22)) ############################################################################### # Section 5 # # The tessellation context class class TessellationContext: """ A context for creating regular tessellations of type {p, q, r} and cusp modulus z = a + b u. # Create a tessellation context for {3,3,6} and z=2+zeta >>> tessContext = TessellationContext(3, 3, 6, 2, 1) """ def __init__(self, p, q, r, a, b): assert (p,q,r) in [(3,3,6), (4,3,6), (5,3,6), (3,4,4)] self.p = p self.q = q self.r = r self.a = a self.b = b self._cachedEdgeCycle = None self._cachedNanotube = None def NewEdgeCycle(self): """ Create a new edge cycle of length 2 * r simplices. This is used a as a building block for nanotubes. >>> tessContext = TessellationContext(3, 3, 6, 2, 1) >>> len(tessContext.NewEdgeCycle()) 12 """ if not self._cachedEdgeCycle: tets = [ Tetrahedron(orientation = (-1)**i) for i in range(2 * self.r) ] for i in range(2 * self.r): GlueFaces(tets[i], tets[(i+1) % (2 * self.r)], 2 + (i % 2)) self._cachedEdgeCycle = set(tets) return deepcopy(self._cachedEdgeCycle) def _NewNanotube(self): tets = set() e = self.NewEdgeCycle() tets |= e tet = FirstPositivelyOrientedTet(e) for m in range(2): otherTet = tet for i in range(self.a): newE = self.NewEdgeCycle() tets |= newE newTet = FirstNegativelyOrientedTet(newE) GlueFaces(otherTet, newTet, 1) newTet.EnforceEdgeCycle(2, 1, 3) otherTet = newTet.GetFaceGluing(3).Rotate(self.r/2, 2, 3) otherTet = otherTet.Rotate(1, 2, 3) for i in range(self.b): newE = self.NewEdgeCycle() tets |= newE newTet = FirstNegativelyOrientedTet(newE) GlueFaces(otherTet, newTet, 1) newTet.EnforceEdgeCycle(2, 1, 3) otherTet = newTet.GetFaceGluing(3).Rotate(self.r/2, 2, 3) otherTet = otherTet.Rotate(-1, 2, 3) Identify(tet, otherTet) tet = tet.Rotate(1, 2, 3) addedTets = FilterCanonical(tets) while addedTets: newTets = set() for tet in addedTets: if tet.GetOrientation() > 0 and not tet.GetFaceGluing(1): newE = self.NewEdgeCycle() newTets |= newE tets |= newE newTet = FirstNegativelyOrientedTet(newE) GlueFaces(tet, newTet, 1) newTet.EnforceEdgeCycle(2, 1, 3) for newTet in newTets: newTet.EnforceEdgeCycle(self.q, 1, 2) addedTets = FilterCanonical(newTets) tets = FilterCanonical(tets) CheckGluingConsistency(tets) for tet in tets: assert tet.Order(1,3) == 2 assert tet.Order(1,2) == self.q assert tet.Order(2,3) == self.r if self.r == 4: assert len(tets) == 8 * ( self.a ** 2 + self.b ** 2) else: assert len(tets) == 12 * ( self.a ** 2 + self.a * self.b + self.b ** 2) return tets def NewNanotube(self): """ Create a new nanotube. >>> tessContext = TessellationContext(3, 3, 6, 2, 1) >>> len(tessContext.NewNanotube()) 84 """ # For speed, nanotubes are cached and a deep copy is created. if not self._cachedNanotube: self._cachedNanotube = self._NewNanotube() sys.setrecursionlimit(10000) return deepcopy(self._cachedNanotube) def EnforceAllEdgeCycles(self, tets): """ For each tet in the given set, enforce that the edges have the given orders. """ while True: keepFixing = False for tet in tets: if tet.EnforceEdgeCycle(2, 0, 2): keepFixing = True if tet.EnforceEdgeCycle(2, 0, 3): keepFixing = True if tet.EnforceEdgeCycle(self.p, 0, 1): keepFixing = True if not keepFixing: return def IsManifold(self, tets): """ Checks whether the triangulation formed by tets is a manifold. >>> tessContext = TessellationContext(3, 3, 6, 2, 0) >>> X = tessContext.PrincipalCongruenceManifold('X') >>> tessContext.IsManifold(X) True >>> Y = tessContext.PrincipalCongruenceManifold('Y') >>> tessContext.IsManifold(Y) False """ orders = (self.p, self.q, self.r, 2, 2, 2) for tet in tets: if not tet.Orders() == orders: return False return True def Extend(self, tets): """ Performs one iteration of step 1 and step 2 of the construction of the universal regular tessellation as described in the paper. It is given the set of tetrahedra that were added during the previous iteration (or a nanotube to get started). It returns the set of tetrahedra that were added in this iteration. >>> tessContext = TessellationContext(3, 3, 6, 3, 0) >>> A = tessContext.NewNanotube() >>> len(A) 108 >>> B = tessContext.Extend(A) >>> len(B) 972 >>> C = tessContext.Extend(B) >>> len(C) 216 >>> D = tessContext.Extend(C) >>> len(D) 0 """ addedTets = set() for i in tets: if i.GetOrientation() > 0 and not i.GetFaceGluing(3).GetFaceGluing(0): newTets = self.NewNanotube() newNanotube = FirstPositivelyOrientedTet(newTets) addedTets |= newTets GlueFaces(i.GetFaceGluing(3), newNanotube, 0) for j in range(2 * self.r): newNanotube.EnforceEdgeCycle(2, 0, 2 + (j % 2)) newNanotube = newNanotube.GetFaceGluing(2 + (j % 2)) self.EnforceAllEdgeCycles(addedTets) return FilterCanonical(addedTets) @staticmethod def GetConnectedNanotube(tet): """ Given a tet, gets all tets forming the nanotube the tet is in. Size of a nanotube >>> tessContext = TessellationContext(3, 3, 6, 2, 0) >>> len(tessContext.NewNanotube()) 48 Nanotube still has same size when constructing manifold >>> X = tessContext.PrincipalCongruenceManifold('X') >>> x = FirstPositivelyOrientedTet(X) >>> len(TessellationContext.GetConnectedNanotube(x)) 48 But in orbifold nanotube has been quotiented by 2 >>> Y = tessContext.PrincipalCongruenceManifold('Y') >>> y = FirstPositivelyOrientedTet(Y) >>> len(TessellationContext.GetConnectedNanotube(y)) 24 """ return GetConnectedComponent(tet, faces = [1, 2, 3]) def UniversalRegularTessellation(self): """ Constructs universal regular tessellation. >>> tessContext = TessellationContext(4, 3, 6, 2, 0) >>> tess = tessContext.UniversalRegularTessellation() >>> len(tess) 768 """ tets = self.NewNanotube() newTets = set([tet for tet in tets]) while newTets: newTets = self.Extend(newTets) tets |= newTets if _verbose: print "Tetrahedra added:", len(newTets) return FilterCanonical(tets) def _IdentityMatrix(self): if self.r == 4: zero = AlgebraicInteger4(0,0) one = AlgebraicInteger4(1,0) else: zero = AlgebraicInteger3(0,0) one = AlgebraicInteger3(1,0) return Matrix2x2(one, zero, zero, one) def _MatrixWhenGoingFaces(self, face1, face2): """ See paper for description """ if self.r == 4: gP = Matrix2x2( AlgebraicInteger4( 0, 0), AlgebraicInteger4( 1, 0), AlgebraicInteger4(-1, 0), AlgebraicInteger4( 1, 0)) gQ = Matrix2x2( AlgebraicInteger4( 0, 1), AlgebraicInteger4( 1, 0), AlgebraicInteger4( 0, 0), AlgebraicInteger4( 1, 0)) gR = Matrix2x2( AlgebraicInteger4( 0, 1), AlgebraicInteger4( 0, 0), AlgebraicInteger4( 0, 0), AlgebraicInteger4( 1, 0)) else: gP = Matrix2x2( AlgebraicInteger3( 0, 0), AlgebraicInteger3( 1, 0), AlgebraicInteger3(-1, 0), AlgebraicInteger3( 1, 0)) gQ = Matrix2x2( AlgebraicInteger3(-1, 1), AlgebraicInteger3( 1, 0), AlgebraicInteger3( 0, 0), AlgebraicInteger3( 1, 0)) gR = Matrix2x2( AlgebraicInteger3( 0, 1), AlgebraicInteger3( 0, 0), AlgebraicInteger3( 0, 0), AlgebraicInteger3( 1, 0)) if face1 == 0: if face2 == 1: return gP if face2 == 2: return gP * gQ if face2 == 3: return gP * gQ * gR if face1 == 1: if face2 == 2: return gQ if face2 == 3: return gQ * gR if face1 == 2: if face2 == 3: return gR def PrincipalCongruenceManifold(self, flavor): """ Constructs a principal congruence manifold or similarly algebraically defined manifold where flavor is "W", "X", "Y". These correspond to the definitions in the paper. >>> tessContext = TessellationContext(3, 3, 6, 1, 0) >>> W = tessContext.PrincipalCongruenceManifold('W') >>> tessContext.IsManifold(W), len(W) (False, 4) >>> X = tessContext.PrincipalCongruenceManifold('X') >>> tessContext.IsManifold(X), len(X) (False, 4) >>> Y = tessContext.PrincipalCongruenceManifold('Y') >>> tessContext.IsManifold(Y), len(Y) (False, 2) >>> tessContext = TessellationContext(3, 3, 6, 1, 1) >>> W = tessContext.PrincipalCongruenceManifold('W') >>> tessContext.IsManifold(W), len(W) (False, 48) >>> X = tessContext.PrincipalCongruenceManifold('X') >>> tessContext.IsManifold(X), len(X) (False, 48) >>> Y = tessContext.PrincipalCongruenceManifold('Y') >>> tessContext.IsManifold(Y), len(Y) (False, 48) >>> tessContext = TessellationContext(3, 3, 6, 2, 0) >>> W = tessContext.PrincipalCongruenceManifold('W') >>> tessContext.IsManifold(W), len(W) (True, 240) >>> X = tessContext.PrincipalCongruenceManifold('X') >>> tessContext.IsManifold(X), len(X) (True, 240) >>> Y = tessContext.PrincipalCongruenceManifold('Y') >>> tessContext.IsManifold(Y), len(Y) (False, 120) >>> tessContext = TessellationContext(3, 3, 6, 2, 1) >>> W = tessContext.PrincipalCongruenceManifold('W') >>> tessContext.IsManifold(W), len(W) (True, 672) >>> X = tessContext.PrincipalCongruenceManifold('X') >>> tessContext.IsManifold(X), len(X) (True, 672) >>> Y = tessContext.PrincipalCongruenceManifold('Y') >>> tessContext.IsManifold(Y), len(Y) (True, 672) >>> tessContext = TessellationContext(3, 3, 6, 2, 2) >>> W = tessContext.PrincipalCongruenceManifold('W') >>> len(W) 2880 >>> X = tessContext.PrincipalCongruenceManifold('X') >>> len(X) 2880 >>> Y = tessContext.PrincipalCongruenceManifold('Y') >>> len(Y) 2880 >>> tessContext = TessellationContext(3, 3, 6, 3, 1) >>> X = tessContext.PrincipalCongruenceManifold('X') >>> len(X) 4368 >>> Y = tessContext.PrincipalCongruenceManifold('Y') >>> len(Y) 2184 >>> tessContext = TessellationContext(3, 3, 6, 4, 0) >>> W = tessContext.PrincipalCongruenceManifold('W') >>> len(W) 7680 >>> X = tessContext.PrincipalCongruenceManifold('X') >>> len(X) 3840 >>> Y = tessContext.PrincipalCongruenceManifold('Y') >>> len(Y) 3840 >>> tessContext = TessellationContext(3, 4, 4, 2, 0) >>> W = tessContext.PrincipalCongruenceManifold('W') >>> tessContext.IsManifold(W), len(W) (True, 192) >>> X = tessContext.PrincipalCongruenceManifold('X') >>> tessContext.IsManifold(X), len(X) (False, 96) >>> Y = tessContext.PrincipalCongruenceManifold('Y') >>> tessContext.IsManifold(Y), len(Y) (False, 96) >>> tessContext = TessellationContext(3, 4, 4, 3, 0) >>> W = tessContext.PrincipalCongruenceManifold('W') >>> len(W) 1440 >>> X = tessContext.PrincipalCongruenceManifold('X') >>> len(X) 1440 >>> Y = tessContext.PrincipalCongruenceManifold('Y') >>> len(Y) 720 """ assert flavor in ['W', 'X', 'Y'] assert self.p == 3, "Not a Bianchi orbifold case" # Switch ring of integers depending on r if self.r == 4: AlgebraicInteger = AlgebraicInteger4 else: AlgebraicInteger = AlgebraicInteger3 # Construct the number z z = AlgebraicInteger(self.a, self.b) # This will be 1 but as type AlgebraicInteger3/4 one = AlgebraicInteger(1, 0) # This will be u u = AlgebraicInteger(0, 1) # Determine the ``ambiguity'' set of algebraic integers. # Two matrices in GL(2,Z[u] / z) are considered equal if they differ # by a constant factor in this set. if flavor == 'W': # In W flavor, this is just +/- 1 ambiguity = set([ one, -one ]) else: # For all other flavors, it is all units in the ring O_d/z ambiguity = UnitsInRingModZ(z) if _verbose: print "Computed units" # Methods for computing the new label for a simplex # For W and X, the label is a pair (matrix, determinant) # where matrix is a matrix in GL(2, Z[u]/z) and determinant # a unit in Z[u]. def NewLabelWX(oldLabel, m): # Multiple the old label with the given matrix and determinant matrix, det = oldLabel matrix, det = matrix * m, det * m.Det() # Bring the label into the form (matrix, 1) or (matrix, u) # by multiplication with a scalar factor while not det in [ one, u ]: matrix = matrix * u det = det * u * u # Return a canonical representative for the equivalence class # of matrices that differ only by a scalar factor in the # ambiguity set return min([ ((matrix * a) % z, det) for a in ambiguity ]) # For Y, the label is just a matrix in GL(2, Z[u]/z) # Computation is similar to W and X labels. def NewLabelY(oldLabel, m): matrix = oldLabel matrix = matrix * m return min([ (matrix * a) % z for a in ambiguity ]) # Pick the appropriate method and label of the root tetrahedron # based on flavor if flavor == 'Y': NewLabel = NewLabelY rootLabel = self._IdentityMatrix() else: NewLabel = NewLabelWX rootLabel = (self._IdentityMatrix(), self._IdentityMatrix().Det()) # A map assigning a label (essentially a matrix) to each simplex labels = {} # A map assigning to each label a simplex inverseLabels = {} def ExtendLabelsFrom(oldTets): """ Given a set of tetrahedra that already have labels, assigning tetrahedra to all connected tetrahedra. Then, identify all tetrahedra which have the same label. """ # List of tetrahedra that need to be processed workingSet = list(oldTets) # A map of a label to a list all tetrahedra that carry this label. # All tetrahedra in each of these lists need to be identified in # the end newInverseLabels = dict([ (k,set([v])) for k, v in inverseLabels.items()]) while workingSet: # Take the next tetrahedron with positive orientation tet = workingSet.pop() if tet.GetOrientation() > 0: # Go along its faces to find other tetrahedron for i1 in range(3): face1 = tet.GetFaceGluing(i1) if face1: for i2 in range(i1 + 1, 4): face2 = face1.GetFaceGluing(i2) # If the other tetrahedron exists and has no # label yet if face2 and not labels.has_key(face2): # Add it to the working set workingSet.append(face2) # Compute the label for the other # tetrahedron m = self._MatrixWhenGoingFaces(i1, i2) newLabel = NewLabel(labels[tet], m) # Assign the label and inverse label labels[face2] = newLabel inverseLabels[newLabel] = face2 # Add it to the list of tets that need to # be identified newInverseLabels[newLabel] = ( set([face2]) | newInverseLabels.get(newLabel, set())) # Identify all tetrahedra which are together in a list. for k, tets in newInverseLabels.items(): first = None for tet in tets: if not first: first = tet else: Identify(first, tet) def ExtendWithLabels(tets): """ Given a set of tetrahedra, attach a new nanotube to each open face, compute the labels for the tetrahedra and identify if necessary. Return the set of newly added tetrahedra. """ addedTets = self.Extend(tets) ExtendLabelsFrom(tets) return addedTets def NanotubeWithLabels(): """ Create the first nanotube, pick a positivley oriented tetrahedron as root and give it with the identity matrix label. Label the other tetrahedra in nanotube as well. """ tets = self.NewNanotube() tet0 = FirstPositivelyOrientedTet(tets) labels[tet0] = rootLabel inverseLabels[rootLabel] = tet0 ExtendLabelsFrom([tet0]) return tets # Get started with a nanotube newTets = NanotubeWithLabels() # All tetrahedra in triangulation tets = FilterCanonical(newTets) while newTets: # Extend by attaching nanotubes to open faces and identify newTets = ExtendWithLabels(newTets) if _verbose: print "Tetrahedra added:", len(newTets) # Keep track of all tetrahedra in triangulation tets = FilterCanonical(tets | newTets) return tets ############################# # Quotient space construction def ApplyWordToTetrahedron(self, tet, word): """ Return the tetrahedron obtained from the given tetrahedron when applying the given word to it. The word is in the finite presentation of the group G in the paper. Each generator in G is encoded as an integer with 1 corresponding to P, 2 to Q, 3 to R. -1 corresponds to P^-1 and so on. New nanotubes are attached automatically if we come along an unglued face when trying to follow the unglued faces. """ otherTet = tet for signedLetter in word: letter = abs(signedLetter) sign = signedLetter / letter self.Extend(TessellationContext.GetConnectedNanotube(otherTet)) otherTet = otherTet.Rotate(sign, letter - 1, letter) return otherTet def IdentifyByWord(self, tet, word): """ Identify the given tetrahedron with its image when applying the word. """ Identify(tet, self.ApplyWordToTetrahedron(tet, word)) def IdentifyByWords(self, tet, words): """ Identify the given tetrahedron with the images when applying words. """ for word in words: self.IdentifyByWord(tet, word) def QuotientByWords(self, words, checkSameCuspModulus = True): """ Construct the quotient space by the action of the subgroup generated by words. words is a list of words encoded as described in ApplyWordToTetrahedron. If checkSameCuspModulus is True, it will check that the map from the universal regular tessellation to the quotient space is a cuspidal covering map. """ # Start with a single nanotube nanotube = self.NewNanotube() nanotubeSize = len(nanotube) # Pick a root tetrahedron, identify it with its images under words tet0 = FirstPositivelyOrientedTet(nanotube) self.IdentifyByWords(tet0, words) # Get the resulting complex we have so far tets = GetConnectedComponent(tet0) newTets = set([tet for tet in tets]) # Extend the resulting complex by adding new nanotubes to open faces # Iterate until no open faces left. while newTets: newTets = self.Extend(newTets) tets |= newTets if _verbose: print "Tetrahedra added:", len(newTets) # The result tets = FilterCanonical(tets) # If asked for, check that each cusp has the same size as a single # nanotube. This will mean that the quotient space has well-defined # cusp modulus and the cusp modulus is equal to the given one. if checkSameCuspModulus: checked = set() for tet in tets: if not tet in checked: nanotube = TessellationContext.GetConnectedNanotube(tet) assert len(nanotube) == nanotubeSize checked |= nanotube self.IsManifold(tets) return tets ############################################################################### # Section 7 # # Special quotients def DodecahedralUniversalQuotientedByC15(): """ Construct a quotient of the universal regular tessellation of type {5, 3, 6} and cusp modulus 2 by the group Z/15. The quotient is tested for having well-defined cusp modulus 2. """ word = [ -1, -1, -2, 1, -2, 1, 1, -3, -1, -1, 3, 2, -1, -1, -2, -3, -3 ] tessContext = TessellationContext(5, 3, 6, 2, 0) return tessContext.QuotientByWords([word]) ############################################################################### # Section 8 # # Test when called from command line def _test(): import doctest doctest.testmod() if __name__ == "__main__": _verbose = False _test()