#***************************************************************************** # Copyright (C) 2007 Carlo Hamalainen , # # Distributed under the terms of the GNU General Public License (GPL) # # This code is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # The full text of the GPL is available at: # # http://www.gnu.org/licenses/ #***************************************************************************** import copy load "dfs_latin.spyx" def back_circulant(n): """ The back-circulant latin square of order n is the Cayley table for (Z_n, +), the integers under addition modulo n. INPUT: n -- int; order of the latin square. EXAMPLE: sage: print back_circulant(5) [0 1 2 3 4] [1 2 3 4 0] [2 3 4 0 1] [3 4 0 1 2] [4 0 1 2 3] """ assert n >= 1 L = matrix(ZZ, n, n) for r in range(n): for c in range(n): L[r, c] = (r + c) % n return L def abelian_2group(s): """ Returns the latin square based on the Cayley table for the abelian 2-group of order 2^s. INPUT: s -- int; order of the latin square will be 2^s. EXAMPLE: sage: print abelian_2group(3) [0 1 2 3 4 5 6 7] [1 0 3 2 5 4 7 6] [2 3 0 1 6 7 4 5] [3 2 1 0 7 6 5 4] [4 5 6 7 0 1 2 3] [5 4 7 6 1 0 3 2] [6 7 4 5 2 3 0 1] [7 6 5 4 3 2 1 0] """ assert s > 0 if s == 1: L = matrix(ZZ, 2, 2) L[0, 0] = 0 L[0, 1] = 1 L[1, 0] = 1 L[1, 1] = 0 return L else: L_prev = abelian_2group(s-1) L = matrix(ZZ, 2^s, 2^s) offset = L.nrows()/2 for r in range(L_prev.nrows()): for c in range(L_prev.ncols()): L[r, c] = L_prev[r, c] L[r+offset, c] = L_prev[r, c] + offset L[r, c+offset] = L_prev[r, c] + offset L[r+offset, c+offset] = L_prev[r, c] return L def empty_latin_square(n): """ Returns an empty latin square of order n. Note that L[r, c] is an empty cell if and only if L[r, c] = -1. INPUT: n -- int; order of the latin square. EXAMPLE: sage: print empty_latin_square(5) [-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] """ assert n >= 1 L = matrix(ZZ, n, n) # fixme: there are so many nested loops like this. Is there a nicer # way to iterate over the cells of a matrix? for r in range(L.nrows()): for c in range(L.ncols()): L[r, c] = -1 return L def top_left_empty_cell(L): """ Returns the least [r, c] such that L[r, c] is an empty cell. INPUT: L -- matrix; a latin square. EXAMPLE: sage: B = back_circulant(5) sage: B[3, 4] = -1 sage: print top_left_empty_cell(B) [3, 4] """ for r in range(L.nrows()): for c in range(L.ncols()): if L[r, c] < 0: return [r, c] return None def permissable_values(L, r, c): """ Find all values that do not appear in row r and column c of the latin square L. If L[r, c] is filled then we return the empty list. INPUT: L -- matrix; a latin square r -- int; row of the latin square c -- int; column of the latin square EXAMPLE: sage: L = back_circulant(5) sage: L[0, 0] = -1 sage: print permissable_values(L, 0, 0) [0] """ if L[r, c] >= 0: return [] assert L.nrows() == L.ncols() n = L.nrows() vals = {} for e in range(n): vals[e] = True for i in range(n): if L[i, c] >= 0: del vals[ L[i, c] ] for j in range(n): if L[r, j] >= 0: try: del vals[ L[r, j] ] except KeyError: # We may have already removed a symbol # in the previous for-loop. pass return vals.keys() def random_empty_cell(P): """ Find an empty cell of P, uniformly at random. INPUT: P -- matrix; a partial latin square RETURNS: [r, c] -- cell such that P[r, c] is empty, or returns None if P is a latin square. """ cells = {} for r in range(P.nrows()): for c in range(P.ncols()): if P[r, c] < 0: cells[ (r,c) ] = True cells = cells.keys() if len(cells) == 0: return None rc = cells[ ZZ.random_element(len(cells)) ] return [rc[0], rc[1]] def dfs_py(P, use_randomised_search, nr_to_find = 0): """ Depth first search to find a completion of the partial latin square P. INPUT: P -- matrix; partial latin square nr_to_find -- int; optional parameter for the maximum number of completions to find. RETURNS: A list of latin squares such that P is contained in each such square. EXAMPLE: sage: P = empty_latin_square(3) sage: P[0, 0] = 0 sage: all_comps = dfs_py(P) sage: some_comps = dfs_py(P, 2) NOTES: There are much quicker ways of finding partial latin square completions (e.g. Knuth's "Dancing Links" algorithm for solving 0-1 matrix exact cover) so this function should only be used for small examples or testing! """ completions = {} def search(P): if use_randomised_search: rc = random_empty_cell(P) else: rc = top_left_empty_cell(P) if rc == None: P.set_immutable() completions[P] = True return r = rc[0] c = rc[1] for e in permissable_values(P, r, c): P2 = copy.deepcopy(P) P2[r, c] = e search(P2) # Make sure that we don't generate too many solutions. if nr_to_find > 0 and len(completions) >= nr_to_find: return search(P) return completions.keys() def is_partial_latin_square(P): """ P is a partial latin square if it is an n by n matrix, and each symbol in [0, 1, ..., n-1] appears at most once in each row, and at most once in each column. """ assert P.nrows() == P.ncols() n = P.nrows() for r in range(n): vals_in_row = {} for c in range(n): e = P[r, c] # Entry out of range 0, 1, ..., n-1: if e >= n: return False # Entry has already appeared in this row: if vals_in_row.has_key(e): return False vals_in_row[e] = True for c in range(n): vals_in_col = {} for r in range(n): e = P[r, c] # Entry out of range 0, 1, ..., n-1: if e >= n: return False # Entry has already appeared in this column: if vals_in_col.has_key(e): return False vals_in_col[e] = True return True def is_latin_square(L): """ L is a partial latin square if it is an n by n matrix, and each symbol in [0, 1, ..., n-1] appears exactly once in each row, and exactly once in each column. """ # We don't allow latin rectangles: if L.nrows() != L.ncols(): return False # Every cell must be filled: if len(filter(lambda x: x >= 0, L.list())) != L.nrows()*L.ncols(): return False # By necessity L must be a partial latin square: if not is_partial_latin_square(L): return False return True def is_uniquely_completable(P): """ Returns True if the partial latin square P has exactly one completion to a latin square. """ return len(dfs_pyrex(P, use_randomised_search = True, nr_to_find = 2)) == 1 def is_completable(P): """ Returns True if the partial latin square can be completed to a latin square. EXAMPLE: The following partial latin square has no completion because there is nowhere that we can place the symbol 0 in the third row: sage: B = empty_latin_square(3) sage: B[0, 0] = 0 sage: B[1, 1] = 0 sage: B[2, 2] = 1 sage: print B sage: print is_completable(B) """ return len(dfs_pyrex(P, use_randomised_search = True, nr_to_find = 1)) > 0 def gcs(L): """ A greedy critical set of a latin square L is found by successively removing elements in a row-wise (bottom-up) manner, checking for unique completion at each step. """ n = L.nrows() G = L.copy() for r in range(n-1, -1, -1): for c in range(n-1, -1, -1): e = G[r, c] G[r, c] = -1 if not is_uniquely_completable(G): G[r, c] = e return G def coin(): """ Simulates a fair coin (returns True or False) using ZZ.random_element(2). """ return ZZ.random_element(2) == 0 def row_containing_sym(L, c, x): """ Given an improper latin square L with L[r1, c] = L[r2, c] = x, return r1 or r2 with equal probability. This is an internal function and should only be used in latin_square_generator(). """ r1 = -1 r2 = -1 for r in range(L.nrows()): if r1 >= 0 and r2 >= 0: break if L[r, c] == x and r1 < 0: r1 = r continue if L[r, c] == x and r2 < 0: r2 = r break assert r1 >= 0 and r2 >= 0 if coin(): return r1 else: return r2 def column_containing_sym(L, r, x): """ Given an improper latin square L with L[r, c1] = L[r, c2] = x, return c1 or c2 with equal probability. This is an internal function and should only be used in latin_square_generator(). """ c1 = -1 c2 = -1 for c in range(L.ncols()): if c1 >= 0 and c2 >= 0: break if L[r, c] == x and c1 < 0: c1 = c continue if L[r, c] == x and c2 < 0: c2 = c break assert c1 >= 0 and c2 >= 0 if coin(): return c1 else: return c2 def next_conjugate(L): """ Permute L[r, c] = e to the conjugate L[c, e] = r We assume that L is an n by n matrix and has values in the range 0, 1, ..., n-1. """ assert L.nrows() == L.ncols() n = L.nrows() C = matrix(ZZ, n, n) for r in range(n): for c in range(n): e = L[r, c] assert e >= 0 and e < n C[c, e] = r return C def latin_square_generator(L_start, check_assertions = False): """ Generator for a sequence of uniformly distributed latin squares, given L_start as the initial latin square. This code implements the Markov chain algorithm of Jacobson and Matthews (1996), see below for the BibTex entry. This generator will never throw the StopIteration exception, so it provides an infinite sequence of latin squares. EXAMPLE: Use the back circulant latin square of order 4 as the initial square and print the next two latin squares given by the Markov chain: sage: g = latin_square_generator(back_circulant(4)) sage: print g.next() sage: print g.next() REFERENCE: @article{MR1410617, AUTHOR = {Jacobson, Mark T. and Matthews, Peter}, TITLE = {Generating uniformly distributed random {L}atin squares}, JOURNAL = {J. Combin. Des.}, FJOURNAL = {Journal of Combinatorial Designs}, VOLUME = {4}, YEAR = {1996}, NUMBER = {6}, PAGES = {405--437}, ISSN = {1063-8539}, MRCLASS = {05B15 (60J10)}, MRNUMBER = {MR1410617 (98b:05021)}, MRREVIEWER = {Lars D{\o}vling Andersen}, } """ if check_assertions: assert is_latin_square(L_start) n = L_start.nrows() r1 = r2 = c1 = c2 = x = y = z = -1 proper = True L = L_start.copy() L_rce = L L_cer = matrix(ZZ, n, n) L_erc = matrix(ZZ, n, n) while True: if proper: if check_assertions: assert is_latin_square(L) ################################# # Update the other two conjugates for r in range(n): for c in range(n): e = L[r, c] L_cer[c, e] = r L_erc[e, r] = c ################################# yield L r1 = ZZ.random_element(n) c1 = ZZ.random_element(n) x = L[r1, c1] y = x while y == x: y = ZZ.random_element(n) # Now find y in row r1 and column c1. if check_assertions: r2 = 0 c2 = 0 while L[r1, c2] != y: c2 += 1 while L[r2, c1] != y: r2 += 1 assert L_erc[y, r1] == c2 assert L_cer[c1, y] == r2 c2 = L_erc[y, r1] r2 = L_cer[c1, y] if check_assertions: assert L[r1, c2] == y if check_assertions: assert L[r2, c1] == y L[r1, c1] = y L[r1, c2] = x L[r2, c1] = x # Now deal with the unknown point. # We want to form z + (y - x) z = L[r2, c2] if z == x: L[r2, c2] = y else: # z and y have positive coefficients # x is the improper term with a negative coefficient proper = False else: # improper square, # L[r2, c2] = y + z - x # y and z are proper while x is the # improper symbol in the cell L[r2, c2]. r1 = row_containing_sym(L, c2, x) c1 = column_containing_sym(L, r2, x) if check_assertions: assert L[r1, c2] == x if check_assertions: assert L[r2, c1] == x # choose one of the proper symbols # uniformly at random (we will use whatever # lands in variable y). if coin(): y, z = z, y # Add/subtract the symbolic difference (y - x) L[r2, c2] = z L[r1, c2] = y L[r2, c1] = y if L[r1, c1] == y: L[r1, c1] = x proper = True else: # got another improper square z = L[r1, c1] x, y = y, x r2 = r1 c2 = c1 # Now we have L[r2, c2] = z+y-x as # usual proper = False # for emphasis def cayley_table(G): """ Construct a latin square on the symbols [0, 1, ..., n-1] for a group with an n by n Cayley table. EXAMPLES: sage: print cayley_table(DihedralGroup(2)) sage: G = gap.Group(PermutationGroupElement((1,2,3))) sage: print cayley_table(G) """ if isinstance(G, sage.interfaces.gap.GapElement): rows = map(lambda x: list(x), list(gap.MultiplicationTable(G))) new_rows = [] for x in rows: new_rows.append(map(lambda x: int(x)-1, x)) return matrix(new_rows) # Otherwise we must have some kind of Sage permutation group object, # such as sage.groups.perm_gps.permgroup.PermutationGroup_generic # or maybe sage.groups.perm_gps.permgroup_named. # We should have a cayley_table() function. try: C = G.cayley_table() except: # We got some other kind of group object? raise NotImplemented L = matrix(ZZ, C.nrows(), C.ncols()) k = 0 entries = {} for r in range(C.nrows()): for c in range(C.ncols()): e_table = C[r, c] if entries.has_key(e_table): L[r, c] = entries[e_table] else: entries[e_table] = k L[r, c] = k k += 1 return L def count_frequencies(n): freq = {} gen = latin_square_generator(back_circulant(n), check_assertions = True) nr_squares = 30000 for _ in range(nr_squares): L = gen.next() M = L.copy() M.set_immutable() try: freq[M] += 1 except KeyError: freq[M] = 1 # fixme: surely there's a built in function to calculate # a frequency distribution in Sage? freq_list = map(lambda x: x/float(nr_squares), freq.values()) return freq_list def alternating_group(m): assert m >= 1 a = tuple(range(1, 2*m+1 + 1)) b = tuple(range(m+1, 0, -1) + range(2*m+2, 3*m+1 + 1)) a = PermutationGroupElement(a) b = PermutationGroupElement(b) c = PermutationGroupElement((a*b)^(-1)) G = PermutationGroup([a, b]) return (a, b, c, G) def pq_group(p, q): assert is_prime(p) assert is_prime(q) assert (q % p) == 1 # beta is a primitive root of the # congruence x^p = 1 mod q F = FiniteField(q) fgen = F.multiplicative_generator() beta = fgen^((q-1)/p) assert beta != 1 assert (beta^p % q) == 1 Q = tuple(range(1, q+1)) P = [] seenValues = {} for i in range(2, q): if seenValues.has_key(i): continue cycle = [] for k in range(p): x = (1 + (i-1)*beta^k) % q if x == 0: x = q seenValues[x] = True cycle.append(x) P.append(tuple(cycle)) G = PermutationGroup([P, Q]) assert G.order() == p*q assert not G.is_abelian() a = PermutationGroupElement(P) b = PermutationGroupElement(Q) c = PermutationGroupElement((a*b)^(-1)) return (a, b, c, PermutationGroup([P, Q])) def p3_group(p): assert is_prime(p) F = gap.new("FreeGroup(3)") a = F.gen(1) b = F.gen(2) c = F.gen(3) rels = [] rels.append( a^p ) rels.append( b^p ) rels.append( c^p ) rels.append( a*b*((b*a*c)^(-1)) ) rels.append( c*a*((a*c)^(-1)) ) rels.append( c*b*((b*c)^(-1)) ) G = F.FactorGroupFpGroupByRels(rels) N = gap.NormalClosure(F, gap.Subgroup(F, rels)) niso = gap.NaturalHomomorphismByNormalSubgroupNC(F, N) a = PermutationGroupElement(gap.Image(niso, a)) b = PermutationGroupElement(gap.Image(niso, b)) c = PermutationGroupElement(gap.Image(niso, c)) G = PermutationGroup([a, b, c]) c = PermutationGroupElement((a*b)^(-1)) return (a, b, c, G) def check_bitrade_properties(a, b, c, G): A = PermutationGroup([a]) B = PermutationGroup([b]) C = PermutationGroup([c]) assert a*b == c^(-1) #assert AB.Size() == AC.Size() == BC.Size() == 1 X = gap.Intersection(gap.Intersection(A, B), C) assert X.Size() == 1 def is_primary_bitrade(a, b, c, G): H = PermutationGroup([a, b, c]) return G.order() == H.order() def coset_representatives(x, G): rcosets = gap.RightCosets(G, PermutationGroup([x])) cosetReps = [] for i in range(1, len(rcosets)+1): cosetReps.append(gap.Representative(rcosets[i])) assert gap.Index(G, PermutationGroup([x])) == len(cosetReps) return cosetReps # The rows of the bitrade are labeled using # the cosets of in G. def bitrade_labels(a, b, c, G): rLabels = cosetRepresentatives(a, G) cLabels = cosetRepresentatives(b, G) eLabels = cosetRepresentatives(c, G) return (rLabels, cLabels, eLabels) # There should be a nicer way to do this, surely? def entry_label(eLabels, C, x): """ Arguments: eLabels --- list of coset representatives of C in G C --- subgroup G of G x --- element of G Returns: integer i such that C*eLabels[i] == C*x """ Cx = gap.RightCoset(C, x) for i in range(len(eLabels)): y = eLabels[i] Cy = gap.RightCoset(C, y) if Cx == Cy: return i assert False def bitrade_pair(a, b, c, G): Tleft = {} Tright = {} (rLabels, cLabels, eLabels) = bitradeLabels(a, b, c, G) T1 = matrix(ZZ, len(rLabels), len(cLabels)) T2 = matrix(ZZ, len(rLabels), len(cLabels)) for i in range(len(rLabels)): for j in range(len(cLabels)): T1[i, j] = -1 T2[i, j] = -1 A = gap.Group([a]) B = gap.Group([b]) C = gap.Group([c]) for i in range(len(rLabels)): for j in range(len(cLabels)): for k in range(len(eLabels)): cosetRow = gap.RightCoset(A, rLabels[i]) cosetCol = gap.RightCoset(B, cLabels[j]) cosetEnt = gap.RightCoset(C, eLabels[k]) #print "---------------------" #print cosetRow #print cosetCol #print cosetEnt g = gap.Intersection(gap.Elements(cosetRow), gap.Elements(cosetCol), gap.Elements(cosetEnt)) if len(gap.Intersection(cosetRow, cosetCol)) == 0: continue if len(gap.Intersection(cosetRow, cosetEnt)) == 0: continue if len(gap.Intersection(cosetCol, cosetEnt)) == 0: continue # At this point we know that a, b, c are the right # generators for a group-based bitrade so we can get # the *unique* element g in the intersection of # cosetRow, cosetCol, and cosetEnt by getting the # intersection of just one pair. g = gap.Intersection(cosetRow, cosetCol) assert len(g) == 1 g = g[1] eIdx = entryLabel(eLabels, C, a^(-1) * g) Tleft[(rLabels[i], cLabels[j])] = eLabels[k] Tright[(rLabels[i], cLabels[j])] = eLabels[eIdx] T1[i, j] = k T2[i, j] = eIdx #return (Tleft, Tright) return (T1, T2) def is_disjoint(Tleft, Tright): assert type(Tleft) == type(Tright) if is_Matrix(Tleft): return is_disjoint_matrix(Tleft, Tright) if isinstance(Tleft, dict): return is_disjoint_dictionary(Tleft, Tright) def is_disjoint_dictionary(Tleft, Tright): #print Tleft.keys() #print Tright.keys() #print Tleft.keys() == Tright.keys() assert Tleft.keys() == Tright.keys() for ij in Tleft.keys(): if Tleft[ij] == Tright[ij]: return False return True def is_disjoint_matrix(Tleft, Tright): for i in range(Tleft.nrows()): for j in range(Tleft.ncols()): if Tleft[i, j] < 0 and Tright[i, j] < 0: continue if Tleft[i, j] == Tright[i, j]: return False return True # examples def example1(): B = back_circulant(6) print B print "Top left entry:", B[0, 0] print "Top-left-most empty cell:", top_left_empty_cell(B) B[3, 4] = -1 print B print "Top-left-most empty cell:", top_left_empty_cell(B) def example2(): A = abelian_2group(3) # abelian 2-group of order 2^3 print A G = gcs(A) print G print "is_uniquely_completable(G) =", is_uniquely_completable(G) print "" G[0, 0] = -1 print G print "is_uniquely_completable(G) =", is_uniquely_completable(G) def example3(): G = gcs(back_circulant(3)) G[0, 0] = -1 print G print "" comps = dfs_pyrex(G, use_randomised_search = True, nr_to_find = 2) for L in comps: print "Completion of G:" print L print "" def example4(): B = back_circulant(4) B[0, 1] = 0 B[1, 0] = -1 print B print "is_partial_latin_square(B) =", is_partial_latin_square(B) def example5(): B = empty_latin_square(3) B[0, 0] = 0 B[1, 1] = 0 B[2, 2] = 1 print B print "is_completable(B) =", is_completable(B) print "" B[2, 2] = 0 print B print "is_completable(B) =", is_completable(B) def example6(): L = abelian_2group(3) print L print "Next conjugate:" print next_conjugate(L) # If we do the conjugate three times we should get back # to the original latin square. print "Conjugate^3(L) == L:", next_conjugate(next_conjugate(next_conjugate(L))) == L def example7(): g = latin_square_generator(back_circulant(4)) print g.next() print "" print g.next() print "" def example8(): P = empty_latin_square(8) # Time some non-randomised searches: for nr_to_find in range(1, 1000, 50): print "aaa", nr_to_find, time dfs_pyrex(P, False, nr_to_find) print "bbb", nr_to_find, time dfs_py(P, False, nr_to_find) # example8() example7()