#!/usr/bin/env python # coding: utf-8 # In[10]: # The file Diagonal.sage contains the routine that is useful to # determine the automorphism group of the 120-cell # This file can be found in the same folder as this one import copy load("Diagonal.sage") maxima_calculus('algebraic: true;') # This helps Sage to deal with fractions phi = (sqrt(5)+1)/2 # If verbose == True, the output contains more information # (used to track what the program is doing, not only to print the result) verbose = False # In[12]: # Some simple functions used to manipulate quaternions # A quaternion here is just a 4-uple of numbers # Quaternion multiplication def quat_mult(z, w): (a,b,c,d,e,f,g,h) = (z[0], z[1], z[2], z[3], w[0], w[1], w[2], w[3]) A = a*e - b*f - c*g - d*h B = a*f + b*e + c*h - d*g C = a*g + c*e + d*f - b*h D = a*h + d*e + b*g - c*f return [A,B,C,D] # Simplify expressions for quaternions (to be used only below - it does not work in general) def simplify_quaternion(q): p = q.copy() for i in range(4): b = q[i] if b != 0 and b != 1/2 and b!= -1/2: c = b.full_simplify() p[i] = c return p # In[14]: # Construct the set S of 6 quaternions in the binary icosahedral group that are closer to 1 S = [[], [], [], [], [], []] S[0] = [phi/2, 0, 1/(2*phi), 1/2] S[1] = [phi/2, 0, 1/(2*phi), -1/2] S[2] = [phi/2, 1/(2*phi), 1/2, 0] S[3] = [phi/2, 1/(2*phi), -1/2, 0] S[4] = [phi/2, 1/2, 0, 1/(2*phi)] S[5] = [phi/2, -1/2, 0, 1/(2*phi)] for i in range(6): S[i] = simplify_quaternion(S[i]) for q in S: if verbose: print (q) # Construct the full binary icosahedral group I of order 120, that is generated by S I = [] for q in S: I.append(q.copy()) sono_tutti = False while sono_tutti == False: for p in I: for q in I: r = quat_mult(p,q) r = simplify_quaternion(r) if I.count(r) == 0: I.append(r) length = len(I) if verbose: print (length, r) if length == 120: sono_tutti = True break if sono_tutti == True: break # The binary icosahedral group I is also the set of the facets of the 120-cell # Sort the elements of I. We automatically get a sorting that respects the # layers 1, 12, 20, 12, 30, 12, 20, 12, 1 of the facets of the 120-cell # The first layer is the south pole [-1, 0, 0, 0], # the second contains the 12 facets adjacent to it, and we end with the north pole [1, 0, 0, 0]. # # Moreover the facets with number i and 119-i are opposite I.sort() for g in I: print(g) # In[15]: # Determine the adjacency matrix G of the 120-cell # G[i][j] = 0 or 1 depending on whether the i-th and j-th vertices of I are far or close G = matrix(120,120) for i in range(120): if verbose: print (i, end = ' ') for j in range(i+1,120): product = sum([I[i][t]*I[j][t] for t in range(4)]) if product == 1/4*sqrt(5) + 1/4: G[i,j] = 1 G[j,i] = 1 # In[7]: # Determine the 14.400 isometries of the 120-cell # This routine is written in Diagonal.sage # It calculates the automorphism of the adjacency graph determined by G # It takes some time isometries = get_isometries(G) # In[9]: # Store the 14.400 isometries on a file import pickle f = open("Isometries", "wb") pickle.dump(isometries, f) # In[16]: # Load the 14.400 isometries from the file # (so we do not need to calculate them every time that we run the program) import pickle g = open("Isometries", "rb") # isometries_copy = pickle.load(g) isometries = pickle.load(g) print(len(isometries)) # In[18]: # Construct all the faces of the 120-cell, as the duals of the flag complex of # the adjacency graph determined by the adjacency matrix G # The facets of the 120-cell are labeled as [0, ..., 119] # and they correspond to the elements in the binary icosahedral group I # faces[i] is the list of all the i-faces of the 120-cell # every i-face f is a (4-i)-uple of numbers in [0, ..., 119] # that correspond to the facets that contain f faces = [[], [], [], []] faces[3] = [i for i in range(120)] for i in range(120): for j in range(i+1,120): if G[i,j] == 1: faces[2].append([i,j]) for i in range(120): for j in range(i+1,120): for k in range(j+1,120): if G[i,j] == 1 and G[i,k] == 1 and G[j,k] == 1: faces[1].append([i,j,k]) for i in range(120): for j in range(i+1,120): for k in range(j+1,120): for l in range(k+1,120): if G[i,j] == 1 and G[i,k] == 1 and G[j,k] == 1 and G[i,l] == 1 and G[j,l] == 1 and G[k,l] == 1: faces[0].append([i,j,k,l]) print (len(faces[0]), len(faces[1]), len(faces[2]), len(faces[3])) for i in range(4): print ("faces of dimension ", i, ": ", faces[i]) # In[21]: # The Davis manifold is obtained from the 120-cell by identifying # opposite facets via translations. # We identify the orbits of the faces after the identification in the Davis manifold # # Warning: we use the fact that the facets i and 119-i are opposite # # The algorithm is quite slow, not very efficient # orbits[i] contains the list of orbits of i-faces orbits = [[], [], [], []] # Orbits of dodecahedral facets for i in range(120): for j in range(i+1, 120): if I[i] == quat_mult(I[j], [-1,0,0,0]): orbits[3].append([i,j]) print ("Orbits of facets: ", orbits[3]) # Orbits of pentagons centers = [] for f in faces[2]: new_center = [(I[f[0]][i] + I[f[1]][i])/2 for i in range(4)] centers.append(new_center) for i in range(len(faces[2])): center = centers[i] adjacent_facets = faces[2][i] gia_trovato = False for o in orbits[2]: for f in o: if set(f) == set(adjacent_facets): gia_trovato = True if gia_trovato == False: new_orbit_centers = [center[:]] keep_going = True if verbose: print (i, center, adjacent_facets) while keep_going: keep_going = False for c in new_orbit_centers: for t in range(2): base_vector = I[adjacent_facets[t]] Fourier = sum([c[i] * base_vector[i] for i in range(4)]) # base_vector is unitary new_c = [c[i] - 2 * Fourier * base_vector[i] for i in range(4)] new_c = simplify_quaternion(new_c) if new_orbit_centers.count(new_c) == 0: new_orbit_centers.append(new_c) keep_going = True if verbose: print (new_c) new_orbit = [] for j in range(len(faces[2])): if new_orbit_centers.count(centers[j]) == 1: new_orbit.append(faces[2][j]) if verbose: print("Orbita:", new_orbit) orbits[2].append(deepcopy(new_orbit)) print ("Orbits of pentagons: ", orbits[2]) # Orbits of edges centers = [] for f in faces[1]: new_center = [(I[f[0]][i] + I[f[1]][i] + I[f[2]][i])/2 for i in range(4)] centers.append(new_center) for i in range(len(faces[1])): center = centers[i] adjacent_facets = faces[1][i] gia_trovato = False for o in orbits[1]: for f in o: if set(f) == set(adjacent_facets): gia_trovato = True if gia_trovato == False: new_orbit_centers = [center[:]] keep_going = True if verbose: print (i, center, adjacent_facets) while keep_going: keep_going = False for c in new_orbit_centers: for t in range(3): base_vector = I[adjacent_facets[t]] Fourier = sum([c[i] * base_vector[i] for i in range(4)]) # base_vector is unitary new_c = [c[i] - 2 * Fourier * base_vector[i] for i in range(4)] new_c = simplify_quaternion(new_c) if new_orbit_centers.count(new_c) == 0: new_orbit_centers.append(new_c) keep_going = True if verbose: print (new_c) new_orbit = [] for j in range(len(faces[1])): if new_orbit_centers.count(centers[j]) == 1: new_orbit.append(faces[1][j]) if verbose: print("Orbita:", new_orbit) orbits[1].append(deepcopy(new_orbit)) print ("Orbits of edges: ", orbits[1]) orbits[0] = [i for i in range(120)] print ("Orbits of vertices: ", orbits[0]) # In[22]: # We get 60 facets, 144 pentagons, 60 edges, and 1 vertex. print("Facets:") print(orbits[3]) print("Pentagons:") print(orbits[2]) print("Edges:") print(orbits[1]) # In[23]: # Store the cellularization on a file import pickle f = open("Orbits", "wb") pickle.dump(orbits, f) # In[24]: # Load the cellularization from the file # (so we do not need to run it every time) import pickle g = open("Orbits", "rb") orbits = pickle.load(g) # In[25]: # For each of the 144 pentagons of the cellularization, write the 5 dodecahedra that contain it # and the 5 edges that are contained in it. # Note that each pentagon has an opposite pentagon with the same data poset = [] # poset[i][0] and poset[i][1] are two 5-uples of numbers in {0, .. , 59} that indicate the 5 dodecahedra # and the 5 edges that contain and are contained in the i-th pentagon. for i in range(144): o = orbits[2][i] p = [[], []] # Determine the dodecahedra that contain the pentagon for j in range(60): trovato = False for t in range(2): for u in range(5): if orbits[2][i][u].count(orbits[3][j][t]) == 1: p[0].append(j) trovato = True break if trovato == True: break # Determine the edges contained in the pentagon for j in range(60): trovato = False for t in range(20): contenuto = True for l in range(2): if orbits[1][j][t].count(orbits[2][i][0][l]) == 0: contenuto = False if contenuto == True: p[1].append(j) trovato = True break poset.append(deepcopy(p)) print (poset) # In[9]: # Recursive function that constructs all the 14400 inside-out isometries. # An inside-out isometry is a permutation p of [0, ..., 59] that indicates a map # that sends the (barycenter of the) i-th dodecahedron to the # (barycenter of the) i-th edge of the cellularization # The dodecahedron i is sent to the edge p[i] # # The isometries are storied in the list isometries[] def inside_out(level = 0, isometry = [], isometries = [], n = 0): if n == 0: n = 60 isometry = [0 for i in range(n)] isometries = [] if level == n: isometries.append(isometry[:]) print (len(isometries), isometry) for i in range(n): if isometry[:level].count(i) == 0: isometry[level] = i if verbose: print (isometry, end='\r', flush = True) is_ok = True for f in range(144): parents = [] for j in range(5): if poset[f][0][j] <= level: parents.append(poset[f][0][j]) trovata = False for g in range(144): contains = True for parent in parents: if poset[g][1].count(isometry[parent]) == 0: contains = False if contains == True: trovata = True break if trovata == False: is_ok = False break if is_ok == True: inside_out(level + 1, isometry, isometries, n) if level == 0: return isometries # In[10]: # Call the recursive function # This produces the list isometries[] of 14400 permutations inside_out() # In[26]: # I write here the first inside-out isometry the routine finds above # Because this is in fact all we need. # Having found this permutation, there is no need to run the routine anymore. IO = [0, 1, 5, 15, 46, 53, 54, 52, 51, 45, 16, 6, 2, 55, 50, 44, 43, 19, 25, 9, 10, 26, 20, 23, 29, 13, 12, 28, 22, 39, 42, 49, 56, 3, 7, 17, 32, 31, 33, 36, 34, 35, 18, 8, 4, 14, 24, 30, 40, 37, 47, 48, 57, 58, 38, 41, 21, 27, 59, 11] print(IO) # In[38]: # Auxiliary routine that checks that the isometry IO is indeed an isometry # This is all we need is_ok = True for p in poset: new_p = [IO[p[0][t]] for t in range(5)] new_p.sort() if verbose: print (p, new_p) quanti = 0 for q in poset: if new_p == q[1]: quanti = quanti + 1 if quanti != 2: is_ok = False print ("IO is not an isometry") break print ("IO is an isometry") # In[40]: # Construct the 72 oriented great circles in the 120-cell as oriented cosets # of the 6 initial order-10 cyclic subgroups # generated by the elements in S. # Each subgroup has 12 cosets # C[12*i+j] is the j-th coset of the i-th subgroup, with i = 0,...,5 and j = 0,...,11 # # Actually, this list could also be deduced more easily from orbits[2], which was written later. C = [] Ccopy = [] # As i goes from 0 to 71, the i-th coset is indicated: # in C[i] as a list of 10 elements of the binary icosahedral group I # in Ccopy[i] as a list of 10 numbers in {0, ..., 119} for i in range(6): # Construct the order-10 cyclic subgroup H generated by q = S[i] q = S[i] p = [1,0,0,0] H = [] for i in range(10): p = quat_mult(p,q) p = simplify_quaternion(p) H.append(p) # Construct the 12 cosets of H CH = [] for p in I: pH = [] for q in H: pq = quat_mult(p,q) pq = simplify_quaternion(pq) pH.append(pq) already_there = False for coset in CH: for a in pH: for b in coset: if a == b: already_there = True if already_there == False: CH.append(pH.copy()) C.append(pH.copy()) if verbose: print (len(C), pH) if len(CH) == 12: break # Reorder cyclically each of the 72 circles so that it starts from a minimmum dodecahedron # Then sort the set of 72 circles for c in C: minimo = 0 for i in range(10): if c[i] < c[minimo]: minimo = i d = [c[(minimo + t) % 10] for t in range(10)] for t in range(10): c[t] = d[t] C.sort() # Translate C into Ccopy, so that the elements of I are denoted as numbers from 0 a 119 for c in C: cnew = [I.index(c[j]) for j in range(10)] Ccopy.append(deepcopy(cnew)) if verbose: print (cnew) for i in range(72): print(Ccopy[i]) # In[47]: # Store Ccopy on a file import pickle f = open("Ccopy", "wb") pickle.dump(Ccopy, f) # In[48]: # Load Ccopy from the file # (so we do not need to calculate it again) import pickle g = open("Ccopy", "rb") Ccopy = pickle.load(g) # In[51]: # Every oriented great circle gives rise to a surface a0 , ..., a71 # (in the paper they are denoted as a1, ..., a72) # # We determine the intersection form Q in the basis (over Q) given by these 72 surfaces. # Q is a 72 x 72 matrix with entries 1, 0, or -1. Given two great circles, we get: # 0 if they intersect, # ±1 if they are disjoint, and the sign is their algebraic intersection number. # We start with a 72x72 matrix M of zeroes and then we modify it M = [([0]*72) for i in range(72)] for i in range(72): for j in range(72): a, b = C[i], C[j] A = matrix([a[0], a[1], b[0], b[1]]) det = A.determinant() M[i][j] = det if M[i][j] != 0: if M[i][j] > 0: M[i][j] = 1 if M[i][j] < 0: M[i][j] = -1 if verbose: print (i, end = ' ') # We copy M to the integer matrix Q Q = Matrix(ZZ, M) # We check that det(Q) = 2^32 d = Q.determinant() if d != 2^72: print ("The determinant is not correct!") # In[52]: # Store the intersection form Q on a file import pickle f = open("Q", "wb") pickle.dump(Q, f) # In[53]: # Load the intersection form Q from the file import pickle g = open("Q", "rb") Q = pickle.load(g) # In[55]: # Check visually that the determinant is 2^72 Q.determinant(), 2^72 # In[56]: # Determine the action of the 14.400 isometries of the 120-cell # on the rational second homology group, with basis a0, ... , a71 that correspond to the oriented great circles. # Each isometry sends each aj to some ±ak. So it is a permutation matrix with arbitrary signs ±1. # # Cisom[i] = [[7,-1], [98,1], [10,-1], ... ] # indicates that the i-th isometries sends a0 to -a7, a98 to a1, etc. # # It turns out that we get only 7.200 different actions. # To these actions we add 7.200 more actions obtained from these by reversing all signs, that is by # composing with -I. This additional action is not induced by an isometry, # but it is still a symmetry of our inequalities. # So at the end we get 14.400 symmetries that preserve the inequalities (hence their solutions set) # We start by listing the 14.400 actions of the isometries Cisom = [] counter = 0 for p in isometries: new_isom = [] for i in range(72): c = Ccopy[i] cimage = [p[c[j]] for j in range(10)] for j in range(72): they_are_equal = True for k in range(10): if Ccopy[j].count(cimage[k]) == 0: they_are_equal = False break if they_are_equal == True: cambio_orientazione = 1 for k in range(10): if cimage[k] == Ccopy[j][0] and cimage[(k+1) % 10] == Ccopy[j][9]: cambio_orientazione = -1 new_isom.append([j,cambio_orientazione]) break if verbose: print (counter, end='\r', flush=True) counter = counter + 1 Cisom.append(deepcopy(new_isom)) # We sort Cisom. # There are duplicates. Change the duplicates by inverting signs. # So we get 14.400 isomorphisms overall. Cisom.sort() for i in range(14399,0,-2): for j in range(72): Cisom[i][j][1] = - Cisom[i][j][1] if verbose: print(len(Cisom)) # Add to each isomorphism p in Cisom an initial list of flags # the flags are the i such that if p[j][0] == i then p[k][0] < i for all k i: is_good = False break if p[j][0] == i: break if is_good: flags.append(i) new_p = [flags[:], p] CCisom.append(deepcopy(new_p)) # Finally, construct the list Isom_efficient, which will be used by the algorithm # Isom_efficient[i] contains the list of all symmetries in CCisom whose initial flags contains i # # These are the symmetries such that all the numbers in the permutation that appear before i are smaller than i # # We avoid constructing Isom_efficient[71] because it would contain all the 144000 isometries # It takes a lot of time to construct it and it would be useless anyway # Hence Isom_efficient[71] is just an empty list by assumption Isom_efficient = [] for i in range(72): new_list = [] if i != 71: # This condition ensures that Isom_efficient[71] is empty (see above) for p in CCisom: if p[0].count(i) == 1: new_list.append(deepcopy(p[1])) Isom_efficient.append(deepcopy(new_list)) print (i, end = ' ') # Removes the identity symmetry everywhere because it is useless for our purposes for i in range(71): Isom_efficient[i].remove([[0, 1], [1, 1], [2, 1], [3, 1], [4, 1], [5, 1], [6, 1], [7, 1], [8, 1], [9, 1], [10, 1], [11, 1], [12, 1], [13, 1], [14, 1], [15, 1], [16, 1], [17, 1], [18, 1], [19, 1], [20, 1], [21, 1], [22, 1], [23, 1], [24, 1], [25, 1], [26, 1], [27, 1], [28, 1], [29, 1], [30, 1], [31, 1], [32, 1], [33, 1], [34, 1], [35, 1], [36, 1], [37, 1], [38, 1], [39, 1], [40, 1], [41, 1], [42, 1], [43, 1], [44, 1], [45, 1], [46, 1], [47, 1], [48, 1], [49, 1], [50, 1], [51, 1], [52, 1], [53, 1], [54, 1], [55, 1], [56, 1], [57, 1], [58, 1], [59, 1], [60, 1], [61, 1], [62, 1], [63, 1], [64, 1], [65, 1], [66, 1], [67, 1], [68, 1], [69, 1], [70, 1], [71, 1]]) # In[58]: # The list Isom_efficient[i] is used to speed up the program, but as i increases the gain in time is smaller # and the time consumed to use it increases, so some balance is needed here. # It is convenient to fix a Cutoff that is determined experimentally, and delete # all the elements in Isom_efficient[i] with i > Cutoff, # because these elements slow down the program more than they speed it up Cutoff = 40 tot = 0 for i in range(Cutoff): tot = tot + len(Isom_efficient[i]) if verbose: print (tot) for i in range(Cutoff, 72): Isom_efficient[i] = [] # In[59]: # Store Isom_efficient on a file import pickle f = open("Isom_efficient", "wb") pickle.dump(Isom_efficient, f) # In[60]: # Load Isom_efficient from the file import pickle g = open("Isom_efficient", "rb") # Isom_efficient_copy = pickle.load(g) Isom_efficient = pickle.load(g) # In[63]: # We determine the action R of the inside-out isometry on the rational homology, with respect to the basis ai. # R is a 72 x 72 matrix # Let OI be the inverse of IO OI = [0 for i in range(60)] for t in range(60): OI[IO[t]] = t print("OI =", OI) # As explained in the paper, each oriented great circle gives rise to two surfaces ai and Ai such that # Q(ai, Aj) = 2 delta(i,j) # The isometry IO sends {ai} to {Aj} without preserving indices, that is it sends ai to Aj where j=P(i) for # some permutation P that we want to determine now poset_semplificato = [] for i in range(144): new_poset = [[], []] for j in range(5): new_poset[0].append(orbits[3][poset[i][0][j]][0]) new_poset[0].append(orbits[3][poset[i][0][j]][1]) for j in range(5): new_poset[1].append(orbits[3][OI[poset[i][1][j]]][0]) new_poset[1].append(orbits[3][OI[poset[i][1][j]]][1]) poset_semplificato.append(deepcopy(new_poset)) P = [0 for i in range(72)] for i in range(72): c = Ccopy[i] for p in poset_semplificato: if set(p[0]) == set(c): for j in range(72): if set(p[1]) == set(Ccopy[j]): P[i] = j print("Action =", P) # We can finally use P to determine the matrix R Qinv = 2 * Q.inverse() R = Matrix(QQ, 72) for i in range(72): for j in range(72): R[i,j] = Qinv[i,P[j]] # The matrix R is actually not quite correct because we have ignored signs. # We can fix them now a posteriori. for k in range(72): for i in range(72): x = R.column(k) * Q * R.column(i) if verbose: print(x, Q[k,i]) if x != Q[k,i]: for j in range(72): R[j,i] = - R[j,i] # Finally we check that R is indeed a Q-isometry. if R.transpose() * Q * R != Q: print ("R is not a Q-isometry") # In[88]: # We now determine the 360 + 360 inequalities produced by the 360 + 360 surfaces b_i and B_i # Each of the 60 dodecahedra in the Davis manifold # contains 6 large hexagonal "equatorial paths" that cut the dodecahedron in two equal pieces. # Each large hexagonal path represents a surface b_i, so these are 60 x 6 = 360 as stated # Each b_i intersects only the 6 elements a_i that go through that dodecahedron, with a sign. # We construct the 360 x 72 matrix B whose element B[i][j] is the intersection between b_i and a_j B = [] for g in I: if g > [0,0,0,0]: # We select only half of the facets, since opposite facets are identified circles = [] # This list will contain the indices of the 6 great circles that cross g for i in range(72): if C[i].count(g) == 1: circles.append(i) for j in range(6): # Each of the 6 circles determines a b_k c = C[circles[j]] # This is the circle that determines b_k cseq = c[(c.index(g) + 1) % 10] v = [0 for i in range(72)] # we want v[i] = for i in range(6): d = C[circles[i]] dseq = d[(d.index(g) + 1) % 10] c_delta = [cseq[i] - g[i] for i in range(4)] d_delta = [dseq[i] - g[i] for i in range(4)] product = sum([c_delta[i] * d_delta[i] for i in range(4)]) if verbose: print (cseq, dseq, product) if product > 0: v[circles[i]] = 1 if product < 0: v[circles[i]] = -1 B.append(v.copy()) # We now construct BB and Beff from the 360 x 72 matrix B. # These are just variations that contain the same data, but expressed in a different way that # will be useful to speed up the routine # # We will then use the more efficient Beff insteand of B # Each line B[i] of B contains precisely 6 non-zero numbers, that are ±1. # The line determines an inequality that involves only 6 variables corresponding to # the positions of the 6 non-zero elements of the line. # # We define BB[i] = [a, B[i]] # where a is the maximum index such that B[i][a] != 0 # Hence the flag a indicates the highest index of the variable involved in the inequality # # Then we reorder the 360 elements in BB, so that the inequalities with smallest flag a are checked first. # # For instance: # B[0] = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0] # BB[0] = [5, [-1, -1, -1, 1, -1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]] BB = [] for b in B: a = 0 for i in range(72): if b[i] != 0: a = i new_b = [a, b[:]] BB.append(new_b) BB.sort() # Beff[a] contains all the lines of B with flag a, for a = 1, ..., 72 # These lines are rewritten using less data following this example: # BB[0] = [5, [-1, -1, -1, 1, -1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]] # becomes # Beff[5][0] = [[0, -1], [1, -1], [2, -1], [3, 1], [4, -1], [5, 1]] Beff = [[] for i in range(72)] for b in BB: new_vector = [] for j in range(72): if b[1][j] != 0: new_vector.append([j, b[1][j]]) Beff[b[0]].append(deepcopy(new_vector)) # Find the 360 inequalities produced by the dual surfaces B_i # These inequalities are obtained from those of B via the action of the inside-out map # We store them as 360 vectors with 72 components in BBdual Bdual = [] for i in range(len(BB)): v = vector(BB[i][1]) w = R.transpose() * v if verbose: print(w) Bdual.append(w) # For instance: # Bdual[0] = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 1, -1, 0, 0, -1, 1, 0, 0, 0, 0, 0, 1, 1, 0, -1, 1, 0, 0, 1, -1, -1, 0, -1, 0, 0, 0, 0, -1, -1, 2, 0) # Since there are more non-zero numbers here, it seems not so useful to try to rewrite this data here # as we did for B. So we keep it as a matrix. # In[89]: # Store Beff and Bdual on a file import pickle f = open("Beff", "wb") pickle.dump(Beff, f) f = open("Bdual", "wb") pickle.dump(Bdual, f) # In[90]: # Load Beff and BDual from the file # So we do not need to calculate them again import pickle g = open("Beff", "rb") Beff = pickle.load(g) g = open("Bdual", "rb") Bdual = pickle.load(g) # In[82]: # # This is the main recursive function that searches for even integral classes x = x_ia_i such that: # # || <= 2, that is |x_i| <= 1 (72 inequalities) # || <= 2, which is calculated using Q (72 inequalities) # || <= 2, which is calculated using Beff (360 inequalities) # || <= 2, which is calculated using Bdual (360 inequalities) # # >= M, this is only auxiliary if one wants to get a shorter output list # # The hypothesis that x is even is checked by requiring that # x_i = -1,0,1 and are all even. # # We set M = 0, but we could put some other value here if needed M = 0 # The counter is the vector vec that encodes the variables x0, ... x71 # vec is an element of {-1,0,1}^72 # We run it lexicographically, but we will cut many dead branches # That is we first set vec = [-1, -10, -10, ..., -10] and ask: is it ok? # (here -10 means that no number -1,0,1 has been chosen yet) # If it is ok, we modify the next number and set vec = [-1, -1, -10, ..., -10] and go on # Then if for instance at some point vec = [-1,-1,-1,-1,-10, ...] is not ok, # we try vec = [-1,-1,-1,0,-10, ...], and go on. # The notion of "being ok" is explained below # level is the number in {0, ..., 71} that indicates that the routine # is going to modify the element vec[level] # n is just an initial flag that is 0 at the beginning and 72 after the first run # answer is the list of all the integral classes found by the algorithm. So it's empty at start. # The inequalities in Beff are already encoded efficiently, but those in Q and Bdual are not, so we now # deal with these. # # To use the inequality in Q and Bdual efficiently we use two matrices called X and missing. # # X and missing are two (72+360) x 72 matrices # # The element X[i] and missing[i] corresponds to the i-th inequality among the 72+360 produced by Q (72) and Bdual (360) # The element X[i][j] encodes the temporary sum of this inequality if one considers only variables <=j # The element missing[i][j] >= 0 encodes the maximum range of variation that this inequality can still have with # the variables >j # # The point is that if abs(X[i][level]) - missing[i][level] > 2 then the inequality |< x, ...>| <= 2 # will never be fulfilled, no matter how we complete the vector vec at the subsequent levels. # Therefore the chosen value for vec[level] is "not ok" and we cut the entire branch. def search_classes(level = 0, vec = [], n = 0, answer = [], X = [], missing = []): if n == 0: # This is the initial part of the recursive algorithm # We set all the variables that will be used to their initial values n = 72 answer = [] vec = [-10 for i in range(72)] X = [[[] for i in range(72)] for j in range(72+360)] missing = [[[] for i in range(72)] for j in range(72+360)] for i in range(72): for j in range(72): temp = [abs(Q[i][k]) for k in range(j+1,72)] missing[i][j] = temp.count(1) for i in range(360): for j in range(72): temp = [abs(Bdual[i][k]) for k in range(j+1,72)] missing[72+i][j] = sum(temp) # This is the general part of the recursive algorithm, that is run many times if level < n: for t in range(-1,2,1): # for t in range(1,-2,-1): # To be used if you want to run the counter in the opposite direction print (vec, end='\r', flush=True) # We set a new value for vec[level] vec[level] = t # We change X[i][level] accordingly for i in range(72): if level == 0: X[i][0] = Q[i][0] * vec[0] else: X[i][level] = X[i][level-1] + Q[i][level] * vec[level] for i in range(360): if level == 0: X[72+i][0] = Bdual[i][0] * vec[0] else: X[72+i][level] = X[72+i][level-1] + Bdual[i][level] * vec[level] is_good = True # Now we put the filters: if something is bad, then is_good becomes False # To be used if you want the counter to start from a 'start' # if level == 20: # start = [-1, -1, -1, -1, -1, -1, -1, -1, 0, -1, 0, -1, 1, 0, -1, 0, 0, -1, 1, 0, 1, 0, 1, 0, -1, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, -1, -1, 0, 0, 1, -1, 1, -1, 0, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -10] # if vec < start: # is_good = False # Check the relevant inequalities in Beff for b in Beff[level]: product = sum([vec[b[i][0]] * b[i][1] for i in range(6)]) if product != -2 and product != 0 and product != 2: is_good = False break if is_good == True: # We now check that the vector vec is minimal under the action of isometries. # Since vec is completed only from 0 to level, the only isometries that make sense # here are those that send {0, ..., level} to itself. # The list Isom_efficient has been constructed to fulfill this requirement: # The list Isom_efficient[level] contains only the isometries of this kind. for p in Isom_efficient[level]: j = 0 while j < 72 and p[j][0] <= level: a = vec[p[j][0]]*p[j][1] if a < vec[j]: is_good = False break if a > vec[j]: break j = j + 1 if is_good == False: break # As explained above we check that abs(X[i][level]) - missing[i][level] > 2 if is_good == True: for i in range(72+360): if abs(X[i][level]) - missing[i][level] > 2: is_good = False break # If is_good is True, we have found no obstructions, and therefore we reiterate the routine if is_good == True: search_classes(level + 1, vec, n, answer, X, missing) if level == 0: return answer if level == n: is_good = True # This is the last step. Here we check again all the filters. w = vector(vec) for i in range(72): product = w*Q[i] if abs(product) > 2: is_good = False return for i in range(360): product = w*Bdual[i] if abs(product) > 2: is_good = False return product = w*Q*w # Here we check that >= M if abs(product) < M: is_good = False return # If there are no obstruction, we print the result and append it to answer if is_good == True: answer.append(vec[:]) print ("found:", abs(product), vec) # In[91]: # This is the main routine L = search_classes()