from collections import Counter import itertools from itertools import groupby #The following takes a list, possibly redundant, and removes redundancies. def unique(a): indices = sorted(range(len(a)), key=a.__getitem__) indices = set(next(it) for k, it in itertools.groupby(indices, key=a.__getitem__)) return [x for i, x in enumerate(a) if i in indices] Id=identity_matrix(4) #Here we make the Gram matrix of the pairing: Q=matrix([[0,2,2,2],[0,0,2,2],[0,0,0,2],[0,0,0,0]]) Q=Q+Q.transpose() #Here we make a list of the generators of the group, together with inverses: s=matrix([[1,2,6,4],[0,-1,-2,-2],[0,2,3,2],[0,0,0,1]]) #translation by a difference of multisections hyp=matrix([[1,0,2,2],[0,1,2,2],[0,0,-1,0],[0,0,0,-1]]) #hyperelliptic involution S=SymmetricGroup(4) S=[x.matrix() for x in S] G1=[x*s*x.transpose() for x in S] G2=[x*hyp*x.transpose() for x in S] G2.extend(G1) H1=[g^-1 for g in G1] H2=[g^-1 for g in G2] G1.extend(H1) G2.extend(H2) GG1=copy(G1) GG1.append(Id) GG2=copy(G2) GG2.append(Id) G2=[x*y for x in GG2 for y in GG2] #Now we are going to build some of the walls for a Dirichlet domain. We start with a vector v=(1,1,1,1), compute part of its orbit g*v0, and then #consider half-spaces containing v and bounded by the hyperplanes orthogonal to g*v-v, where 'orthogonal' means with respect to Q. v=matrix([[1],[1],[1],[1]]) HYP1=[(g*v -v).transpose()*Q for g in G1] HYP2=[(g*v -v).transpose()*Q for g in G2] Z=matrix([0]) R1=[Z.augment(h) for h in HYP1] R2=[Z.augment(h) for h in HYP2] R1=[list(R1[i])[0] for i in range(len(R1))] R2=[list(R2[i])[0] for i in range(len(R2))] P1=Polyhedron(ieqs=R1) P2=Polyhedron(ieqs=R2) WALLS1=P1.Hrepresentation() WALLS2=P2.Hrepresentation() #Output of WALLS2: #(An inequality (0, 0, 1, 1) x + 0 >= 0, # An inequality (0, 1, 0, 1) x + 0 >= 0, # An inequality (0, 1, 1, 0) x + 0 >= 0, # An inequality (1, 0, 0, 1) x + 0 >= 0, # An inequality (1, 0, 1, 0) x + 0 >= 0, # An inequality (1, 1, 0, 0) x + 0 >= 0) D=Cone([(0,0,1,1),(0,1,0,1),(0,1,1,0),(1,0,0,1),(1,0,1,0),(1,1,0,0)]) C=D.dual() Generators=C.rays() #Outout of Generators: #M( 1, 1, -1, 1), #M( 1, 0, 0, 0), #M( 1, 1, 1, -1), #M( 0, 1, 0, 0), #M(-1, 1, 1, 1), #M( 0, 0, 0, 1), #M( 1, -1, 1, 1), #M( 0, 0, 1, 0) #in 4-d lattice M