Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagelib
Path: blob/master/sage/combinat/matrices/latin.py
4069 views
1
r"""
2
Latin Squares
3
4
A *latin square* of order `n` is an `n \times n` array such that
5
each symbol `s \in \{ 0, 1, \dots, n-1\}` appears precisely once in each
6
row, and precisely once in each column. A *partial latin square* of
7
order `n` is an `n \times n` array such that
8
each symbol `s \in \{ 0, 1, \dots, n-1\}` appears at most once in each
9
row, and at most once in each column. Empty cells are denoted by `-1`.
10
A latin square `L` is a
11
*completion* of a partial latin square `P` if `P \subseteq L`. If
12
`P` completes to just `L` then `P` *has unique completion*.
13
14
A *latin bitrade* `(T_1,\, T_2)` is a pair of partial
15
latin squares such that:
16
17
#. `\{ (i,\,j) \mid (i,\,j,\,k) \in T_1 \mbox{ for some symbol $k$} \} = \{ (i,\,j) \mid (i,\,j,\,k') \in T_2 \mbox{ for some symbol $k'$} \};`
18
19
#. for each `(i,\,j,\,k) \in T_1` and `(i,\,j,\,k') \in T_2`,
20
`k \neq k'`;
21
22
#. the symbols appearing in row `i` of `T_1` are the same as those of
23
row `i` of `T_2`; the symbols appearing in column `j` of `T_1` are
24
the same as those of column `j` of `T_2`.
25
26
27
Intuitively speaking, a bitrade gives the difference between two latin
28
squares, so if `(T_1,\, T_2)` is a bitrade
29
for the pair of latin squares `(L_1,\, L_2)`, then
30
`L1 = (L2 \setminus T_1) \cup T_2`
31
and
32
`L2 = (L1 \setminus T_2) \cup T_1`.
33
34
This file contains
35
36
37
#. LatinSquare class definition;
38
39
#. some named latin squares (back circulant, forward circulant, abelian
40
`2`-group);
41
42
#. functions is\_partial\_latin\_square and is\_latin\_square to test
43
if a LatinSquare object satisfies the definition of a latin square
44
or partial latin square, respectively;
45
46
#. tests for completion and unique completion (these use the C++
47
implementation of Knuth's dancing links algorithm to solve the
48
problem as a instance of `0-1` matrix exact cover);
49
50
#. functions for calculating the `\tau_i` representation of a bitrade
51
and the genus of the associated hypermap embedding;
52
53
#. Markov chain of Jacobson and Matthews (1996) for generating latin
54
squares uniformly at random (provides a generator interface);
55
56
#. a few examples of `\tau_i` representations of bitrades constructed
57
from the action of a group on itself by right multiplication,
58
functions for converting to a pair of LatinSquare objects.
59
60
EXAMPLES::
61
62
sage: from sage.combinat.matrices.latin import *
63
sage: B = back_circulant(5)
64
sage: B
65
[0 1 2 3 4]
66
[1 2 3 4 0]
67
[2 3 4 0 1]
68
[3 4 0 1 2]
69
[4 0 1 2 3]
70
sage: B.is_latin_square()
71
True
72
sage: B[0, 1] = 0
73
sage: B.is_latin_square()
74
False
75
76
sage: (a, b, c, G) = alternating_group_bitrade_generators(1)
77
sage: (T1, T2) = bitrade_from_group(a, b, c, G)
78
sage: T1
79
[ 0 2 -1 1]
80
[-1 0 1 3]
81
[ 3 -1 0 2]
82
[ 1 3 2 -1]
83
sage: T2
84
[ 1 0 -1 2]
85
[-1 3 0 1]
86
[ 0 -1 2 3]
87
[ 3 2 1 -1]
88
sage: T1.nr_filled_cells()
89
12
90
sage: genus(T1, T2)
91
1
92
93
To do:
94
95
96
#. Latin squares with symbols from a ring instead of the integers
97
`\{ 0, 1, \dots, n-1 \}`.
98
99
#. Isotopism testing of latin squares and bitrades via graph
100
isomorphism (nauty?).
101
102
#. Combinatorial constructions for bitrades.
103
104
105
AUTHORS:
106
107
- Carlo Hamalainen (2008-03-23): initial version
108
109
110
TESTS::
111
112
sage: L = elementary_abelian_2group(3)
113
sage: L == loads(dumps(L))
114
True
115
116
"""
117
118
#*****************************************************************************
119
# Copyright (C) 2008 Carlo Hamalainen <[email protected]>,
120
#
121
# Distributed under the terms of the GNU General Public License (GPL)
122
#
123
# This code is distributed in the hope that it will be useful,
124
# but WITHOUT ANY WARRANTY; without even the implied warranty of
125
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
126
# General Public License for more details.
127
#
128
# The full text of the GPL is available at:
129
#
130
# http://www.gnu.org/licenses/
131
#*****************************************************************************
132
133
from sage.matrix.all import matrix
134
from sage.rings.all import ZZ
135
from sage.rings.all import Integer
136
from sage.matrix.matrix_integer_dense import Matrix_integer_dense
137
from sage.groups.perm_gps.permgroup_element import PermutationGroupElement
138
from sage.interfaces.gap import GapElement
139
from sage.combinat.permutation import Permutation
140
from sage.interfaces.gap import gap
141
from sage.groups.perm_gps.permgroup import PermutationGroup
142
from sage.rings.arith import is_prime
143
from sage.rings.finite_rings.constructor import FiniteField
144
from sage.misc.misc import uniq
145
from sage.misc.flatten import flatten
146
147
#load "dancing_links.spyx"
148
#load "dancing_links.sage"
149
150
from dlxcpp import DLXCPP
151
152
class LatinSquare:
153
def __init__(self, *args):
154
"""
155
Latin squares.
156
157
This class implements a latin square of order n with rows and
158
columns indexed by the set 0, 1, ..., n-1 and symbols from the same
159
set. The underlying latin square is a matrix(ZZ, n, n). If L is a
160
latin square, then the cell at row r, column c is empty if and only
161
if L[r, c] < 0. In this way we allow partial latin squares and can
162
speak of completions to latin squares, etc.
163
164
There are two ways to declare a latin square:
165
166
Empty latin square of order n::
167
168
sage: n = 3
169
sage: L = LatinSquare(n)
170
sage: L
171
[-1 -1 -1]
172
[-1 -1 -1]
173
[-1 -1 -1]
174
175
Latin square from a matrix::
176
177
sage: M = matrix(ZZ, [[0, 1], [2, 3]])
178
sage: LatinSquare(M)
179
[0 1]
180
[2 3]
181
"""
182
183
if len(args) == 1 and (type(args[0]) == Integer or type(args[0]) == int):
184
self.square = matrix(ZZ, args[0], args[0])
185
self.clear_cells()
186
elif len(args) == 2 and (type(args[0]) == Integer or type(args[0]) == int) and (type(args[1]) == Integer or type(args[1]) == int):
187
self.square = matrix(ZZ, args[0], args[1])
188
self.clear_cells()
189
elif len(args) == 1 and type(args[0]) == Matrix_integer_dense:
190
self.square = args[0]
191
else:
192
raise NotImplemented
193
194
def dumps(self):
195
"""
196
Since the latin square class doesn't hold any other private
197
variables we just call dumps on self.square:
198
199
EXAMPLES::
200
201
sage: from sage.combinat.matrices.latin import *
202
sage: back_circulant(2) == loads(dumps(back_circulant(2)))
203
True
204
"""
205
206
return dumps(self.square)
207
208
def __str__(self):
209
"""
210
The string representation of a latin square is the same as the
211
underlying matrix.
212
213
EXAMPLES::
214
215
sage: print LatinSquare(matrix(ZZ, [[0, 1], [2, 3]])).__str__()
216
[0 1]
217
[2 3]
218
"""
219
return self.square.__str__()
220
221
def __repr__(self):
222
"""
223
The representation of a latin square is the same as the underlying
224
matrix.
225
226
EXAMPLES::
227
228
sage: print LatinSquare(matrix(ZZ, [[0, 1], [2, 3]])).__repr__()
229
[0 1]
230
[2 3]
231
"""
232
return self.square.__str__()
233
return self.square.__repr__()
234
235
def __getitem__(self, rc):
236
"""
237
If L is a LatinSquare then this method allows us to evaluate L[r,
238
c].
239
240
EXAMPLES::
241
242
sage: from sage.combinat.matrices.latin import *
243
sage: B = back_circulant(3)
244
sage: B[1, 1]
245
2
246
"""
247
248
r = rc[0]
249
c = rc[1]
250
251
return self.square[r, c]
252
253
def __setitem__(self, rc, val):
254
"""
255
If L is a LatinSquare then this method allows us to set L[r, c].
256
257
EXAMPLES::
258
259
sage: from sage.combinat.matrices.latin import *
260
sage: B = back_circulant(3)
261
sage: B[1, 1] = 10
262
sage: B[1, 1]
263
10
264
"""
265
266
r = rc[0]
267
c = rc[1]
268
269
self.square[r, c] = val
270
271
def set_immutable(self):
272
"""
273
A latin square is immutable if the underlying matrix is immutable.
274
275
EXAMPLES::
276
277
sage: L = LatinSquare(matrix(ZZ, [[0, 1], [2, 3]]))
278
sage: L.set_immutable()
279
sage: {L : 0} # this would fail without set_immutable()
280
{[0 1]
281
[2 3]: 0}
282
"""
283
284
self.square.set_immutable()
285
286
def __hash__(self):
287
"""
288
The hash of a latin square is precisely the hash of the underlying
289
matrix.
290
291
EXAMPLES::
292
293
sage: L = LatinSquare(matrix(ZZ, [[0, 1], [2, 3]]))
294
sage: L.set_immutable()
295
sage: L.__hash__()
296
12
297
"""
298
299
return self.square.__hash__()
300
301
def __eq__(self, Q):
302
"""
303
Two latin squares are equal if the underlying matrices are equal.
304
305
EXAMPLES::
306
307
sage: A = LatinSquare(matrix(ZZ, [[0, 1], [2, 3]]))
308
sage: B = LatinSquare(matrix(ZZ, [[0, 4], [2, 3]]))
309
sage: A == B
310
False
311
sage: B[0, 1] = 1
312
sage: A == B
313
True
314
"""
315
316
return self.square == Q.square
317
318
def __copy__(self):
319
"""
320
To copy a latin square we must copy the underlying matrix.
321
322
EXAMPLES::
323
324
sage: A = LatinSquare(matrix(ZZ, [[0, 1], [2, 3]]))
325
sage: B = copy(A)
326
sage: B
327
[0 1]
328
[2 3]
329
"""
330
C = LatinSquare(self.square.nrows(), self.square.ncols())
331
from copy import copy
332
C.square = copy(self.square)
333
return C
334
335
def clear_cells(self):
336
"""
337
Mark every cell in self as being empty.
338
339
EXAMPLES::
340
341
sage: A = LatinSquare(matrix(ZZ, [[0, 1], [2, 3]]))
342
sage: A.clear_cells()
343
sage: A
344
[-1 -1]
345
[-1 -1]
346
"""
347
348
for r in range(self.square.nrows()):
349
for c in range(self.square.ncols()):
350
self.square[r, c] = -1;
351
352
def nrows(self):
353
"""
354
Number of rows in the latin square.
355
356
EXAMPLES::
357
358
sage: LatinSquare(3).nrows()
359
3
360
"""
361
362
return self.square.nrows()
363
364
def ncols(self):
365
"""
366
Number of columns in the latin square.
367
368
EXAMPLES::
369
370
sage: LatinSquare(3).ncols()
371
3
372
"""
373
return self.square.ncols()
374
375
def row(self, x):
376
"""
377
Returns row x of the latin square.
378
379
EXAMPLES::
380
381
sage: from sage.combinat.matrices.latin import *
382
sage: back_circulant(3).row(0)
383
(0, 1, 2)
384
"""
385
386
return self.square.row(x)
387
388
def column(self, x):
389
"""
390
Returns column x of the latin square.
391
392
EXAMPLES::
393
394
sage: from sage.combinat.matrices.latin import *
395
sage: back_circulant(3).column(0)
396
(0, 1, 2)
397
"""
398
return self.square.column(x)
399
400
def list(self):
401
"""
402
Convert the latin square into a list, in a row-wise manner.
403
404
EXAMPLES::
405
406
sage: from sage.combinat.matrices.latin import *
407
sage: back_circulant(3).list()
408
[0, 1, 2, 1, 2, 0, 2, 0, 1]
409
"""
410
411
return self.square.list()
412
413
def nr_filled_cells(self):
414
"""
415
Returns the number of filled cells (i.e. cells with a positive
416
value) in the partial latin square self.
417
418
EXAMPLES::
419
420
sage: from sage.combinat.matrices.latin import *
421
sage: LatinSquare(matrix([[0, -1], [-1, 0]])).nr_filled_cells()
422
2
423
"""
424
425
s = 0
426
for r in range(self.nrows()):
427
for c in range(self.ncols()):
428
if self[r, c] >= 0: s += 1
429
return s
430
431
def actual_row_col_sym_sizes(self):
432
"""
433
Bitrades sometimes end up in partial latin squares with unused
434
rows, columns, or symbols. This function works out the actual
435
number of used rows, columns, and symbols.
436
437
.. warning::
438
439
We assume that the unused rows/columns occur in the lower
440
right of self, and that the used symbols are in the range
441
{0, 1, ..., m} (no holes in that list).
442
443
EXAMPLE::
444
445
sage: from sage.combinat.matrices.latin import *
446
sage: B = back_circulant(3)
447
sage: B[0,2] = B[1,2] = B[2,2] = -1
448
sage: B[0,0] = B[2,1] = -1
449
sage: B
450
[-1 1 -1]
451
[ 1 2 -1]
452
[ 2 -1 -1]
453
sage: B.actual_row_col_sym_sizes()
454
(3, 2, 2)
455
"""
456
457
row_max = self.nrows()
458
col_max = self.ncols()
459
sym_max = self.nr_distinct_symbols()
460
461
while self.is_empty_row(row_max-1): row_max -= 1
462
while self.is_empty_column(col_max-1): col_max -= 1
463
464
return row_max, col_max, sym_max
465
466
def is_empty_column(self, c):
467
"""
468
Checks if column c of the partial latin square self is empty.
469
470
EXAMPLES::
471
472
sage: from sage.combinat.matrices.latin import *
473
sage: L = back_circulant(4)
474
sage: L.is_empty_column(0)
475
False
476
sage: L[0,0] = L[1,0] = L[2,0] = L[3,0] = -1
477
sage: L.is_empty_column(0)
478
True
479
"""
480
481
return uniq(self.column(c)) == [-1]
482
483
def is_empty_row(self, r):
484
"""
485
Checks if row r of the partial latin square self is empty.
486
487
EXAMPLES::
488
489
sage: from sage.combinat.matrices.latin import *
490
sage: L = back_circulant(4)
491
sage: L.is_empty_row(0)
492
False
493
sage: L[0,0] = L[0,1] = L[0,2] = L[0,3] = -1
494
sage: L.is_empty_row(0)
495
True
496
"""
497
498
return uniq(self.row(r)) == [-1]
499
500
def nr_distinct_symbols(self):
501
"""
502
Returns the number of distinct symbols in the partial latin square
503
self.
504
505
EXAMPLE::
506
507
sage: from sage.combinat.matrices.latin import *
508
sage: back_circulant(5).nr_distinct_symbols()
509
5
510
sage: L = LatinSquare(10)
511
sage: L.nr_distinct_symbols()
512
0
513
sage: L[0, 0] = 0
514
sage: L[0, 1] = 1
515
sage: L.nr_distinct_symbols()
516
2
517
"""
518
519
symbols = uniq(flatten(map(lambda x: list(x), list(self.square))))
520
symbols = filter(lambda x: x >= 0, symbols)
521
522
return len(symbols)
523
524
def apply_isotopism(self, row_perm, col_perm, sym_perm):
525
"""
526
An isotopism is a permutation of the rows, columns, and symbols of
527
a partial latin square self. Use isotopism() to convert a tuple
528
(indexed from 0) to a Permutation object.
529
530
EXAMPLES::
531
532
sage: from sage.combinat.matrices.latin import *
533
sage: B = back_circulant(5)
534
sage: B
535
[0 1 2 3 4]
536
[1 2 3 4 0]
537
[2 3 4 0 1]
538
[3 4 0 1 2]
539
[4 0 1 2 3]
540
sage: alpha = isotopism((0,1,2,3,4))
541
sage: beta = isotopism((1,0,2,3,4))
542
sage: gamma = isotopism((2,1,0,3,4))
543
sage: B.apply_isotopism(alpha, beta, gamma)
544
[3 4 2 0 1]
545
[0 2 3 1 4]
546
[1 3 0 4 2]
547
[4 0 1 2 3]
548
[2 1 4 3 0]
549
"""
550
551
#Q = matrix(ZZ, self.nrows(), self.ncols())
552
Q = LatinSquare(self.nrows(), self.ncols())
553
554
for r in range(self.nrows()):
555
for c in range(self.ncols()):
556
try:
557
if self[r, c] < 0: s2 = -1
558
else: s2 = sym_perm[self[r, c]] - 1
559
except IndexError:
560
s2 = self[r, c] # we must be leaving the symbol fixed?
561
562
Q[row_perm[r]-1, col_perm[c]-1] = s2
563
564
return Q
565
566
def filled_cells_map(self):
567
"""
568
Number the filled cells of self with integers from {1, 2, 3, ...}
569
570
INPUT:
571
572
573
- ``self`` - Partial latin square self (empty cells
574
have negative values)
575
576
577
OUTPUT: A dictionary cells_map where cells_map[(i,j)] = m means
578
that (i,j) is the m-th filled cell in P, while cells_map[m] =
579
(i,j).
580
581
EXAMPLES::
582
583
sage: from sage.combinat.matrices.latin import *
584
sage: (a, b, c, G) = alternating_group_bitrade_generators(1)
585
sage: (T1, T2) = bitrade_from_group(a, b, c, G)
586
sage: T1.filled_cells_map()
587
{1: (0, 0), 2: (0, 1), 3: (0, 3), 4: (1, 1), 5: (1, 2), 6: (1, 3), 7: (2, 0), 8: (2, 2), 9: (2, 3), 10: (3, 0), 11: (3, 1), 12: (3, 2), (1, 3): 6, (0, 3): 3, (1, 2): 5, (3, 0): 10, (2, 2): 8, (1, 1): 4, (3, 2): 12, (0, 0): 1, (2, 3): 9, (0, 1): 2, (3, 1): 11, (2, 0): 7}
588
"""
589
590
cells_map = {}
591
k = 1
592
593
for r in range(self.nrows()):
594
for c in range(self.ncols()):
595
e = self[r, c]
596
597
if e < 0: continue
598
599
cells_map[ (r,c) ] = k
600
cells_map[k] = (r,c)
601
602
k += 1
603
604
return cells_map
605
606
def top_left_empty_cell(self):
607
"""
608
Returns the least [r, c] such that self[r, c] is an empty cell. If
609
all cells are filled then we return None.
610
611
INPUT:
612
613
614
- ``self`` - LatinSquare
615
616
617
EXAMPLES::
618
619
sage: from sage.combinat.matrices.latin import *
620
sage: B = back_circulant(5)
621
sage: B[3, 4] = -1
622
sage: B.top_left_empty_cell()
623
[3, 4]
624
"""
625
626
for r in range(self.nrows()):
627
for c in range(self.ncols()):
628
if self[r, c] < 0:
629
return [r, c]
630
631
return None
632
633
def is_partial_latin_square(self):
634
"""
635
self is a partial latin square if it is an n by n matrix, and each
636
symbol in [0, 1, ..., n-1] appears at most once in each row, and at
637
most once in each column.
638
639
EXAMPLES::
640
641
sage: from sage.combinat.matrices.latin import *
642
sage: LatinSquare(4).is_partial_latin_square()
643
True
644
sage: back_circulant(3).gcs().is_partial_latin_square()
645
True
646
sage: back_circulant(6).is_partial_latin_square()
647
True
648
"""
649
650
assert self.nrows() == self.ncols()
651
652
n = self.nrows()
653
654
for r in range(n):
655
vals_in_row = {}
656
657
for c in range(n):
658
e = self[r, c]
659
660
if e < 0: continue
661
662
# Entry out of range 0, 1, ..., n-1:
663
if e >= n: return False
664
665
# Entry has already appeared in this row:
666
if vals_in_row.has_key(e): return False
667
668
vals_in_row[e] = True
669
670
for c in range(n):
671
vals_in_col = {}
672
673
for r in range(n):
674
e = self[r, c]
675
676
if e < 0: continue
677
678
# Entry out of range 0, 1, ..., n-1:
679
if e >= n: return False
680
681
# Entry has already appeared in this column:
682
if vals_in_col.has_key(e): return False
683
684
vals_in_col[e] = True
685
686
return True
687
688
def is_latin_square(self):
689
"""
690
self is a latin square if it is an n by n matrix, and each symbol
691
in [0, 1, ..., n-1] appears exactly once in each row, and exactly
692
once in each column.
693
694
EXAMPLES::
695
696
sage: from sage.combinat.matrices.latin import *
697
sage: elementary_abelian_2group(4).is_latin_square()
698
True
699
700
::
701
702
sage: forward_circulant(7).is_latin_square()
703
True
704
"""
705
706
# We don't allow latin rectangles:
707
if self.nrows() != self.ncols():
708
return False
709
710
# Every cell must be filled:
711
if len(filter(lambda x: x >= 0, self.list())) != self.nrows()*self.ncols():
712
return False
713
714
# By necessity self must be a partial latin square:
715
if not self.is_partial_latin_square():
716
return False
717
718
return True
719
720
def permissable_values(self, r, c):
721
"""
722
Find all values that do not appear in row r and column c of the
723
latin square self. If self[r, c] is filled then we return the empty
724
list.
725
726
INPUT:
727
728
729
- ``self`` - LatinSquare
730
731
- ``r`` - int; row of the latin square
732
733
- ``c`` - int; column of the latin square
734
735
736
EXAMPLES::
737
738
sage: from sage.combinat.matrices.latin import *
739
sage: L = back_circulant(5)
740
sage: L[0, 0] = -1
741
sage: L.permissable_values(0, 0)
742
[0]
743
"""
744
745
if self[r, c] >= 0:
746
return []
747
748
assert self.nrows() == self.ncols()
749
750
n = self.nrows()
751
752
vals = {}
753
for e in range(n):
754
vals[e] = True
755
756
for i in range(n):
757
if self[i, c] >= 0:
758
del vals[ self[i, c] ]
759
760
for j in range(n):
761
if self[r, j] >= 0:
762
try:
763
del vals[ self[r, j] ]
764
except KeyError:
765
# We may have already removed a symbol
766
# in the previous for-loop.
767
pass
768
769
return vals.keys()
770
771
def random_empty_cell(self):
772
"""
773
Find an empty cell of self, uniformly at random.
774
775
INPUT:
776
777
- ``self`` - LatinSquare
778
779
OUTPUT:
780
781
- ``[r, c]`` - cell such that self[r, c] is empty, or returns
782
None if self is a (full) latin square.
783
784
EXAMPLES::
785
786
sage: from sage.combinat.matrices.latin import *
787
sage: P = back_circulant(2)
788
sage: P[1,1] = -1
789
sage: P.random_empty_cell()
790
[1, 1]
791
"""
792
793
cells = {}
794
795
for r in range(self.nrows()):
796
for c in range(self.ncols()):
797
if self[r, c] < 0:
798
cells[ (r,c) ] = True
799
800
cells = cells.keys()
801
802
if len(cells) == 0: return None
803
804
rc = cells[ ZZ.random_element(len(cells)) ]
805
806
return [rc[0], rc[1]]
807
808
def is_uniquely_completable(self):
809
"""
810
Returns True if the partial latin square self has exactly one
811
completion to a latin square. This is just a wrapper for the
812
current best-known algorithm, Dancing Links by Knuth. See
813
dancing_links.spyx
814
815
EXAMPLES::
816
817
sage: from sage.combinat.matrices.latin import *
818
sage: back_circulant(4).gcs().is_uniquely_completable()
819
True
820
821
::
822
823
sage: G = elementary_abelian_2group(3).gcs()
824
sage: G.is_uniquely_completable()
825
True
826
827
::
828
829
sage: G[0, 0] = -1
830
sage: G.is_uniquely_completable()
831
False
832
"""
833
834
return self.dlxcpp_has_unique_completion()
835
836
def is_completable(self):
837
"""
838
Returns True if the partial latin square can be completed to a
839
latin square.
840
841
EXAMPLES:
842
843
The following partial latin square has no completion because there
844
is nowhere that we can place the symbol 0 in the third row::
845
846
sage: B = LatinSquare(3)
847
848
::
849
850
sage: B[0, 0] = 0
851
sage: B[1, 1] = 0
852
sage: B[2, 2] = 1
853
854
::
855
856
sage: B
857
[ 0 -1 -1]
858
[-1 0 -1]
859
[-1 -1 1]
860
861
::
862
863
sage: B.is_completable()
864
False
865
866
::
867
868
sage: B[2, 2] = 0
869
sage: B.is_completable()
870
True
871
"""
872
873
return len(dlxcpp_find_completions(self, nr_to_find = 1)) > 0
874
875
def gcs(self):
876
"""
877
A greedy critical set of a latin square self is found by
878
successively removing elements in a row-wise (bottom-up) manner,
879
checking for unique completion at each step.
880
881
EXAMPLES::
882
883
sage: from sage.combinat.matrices.latin import *
884
sage: A = elementary_abelian_2group(3)
885
sage: G = A.gcs()
886
sage: A
887
[0 1 2 3 4 5 6 7]
888
[1 0 3 2 5 4 7 6]
889
[2 3 0 1 6 7 4 5]
890
[3 2 1 0 7 6 5 4]
891
[4 5 6 7 0 1 2 3]
892
[5 4 7 6 1 0 3 2]
893
[6 7 4 5 2 3 0 1]
894
[7 6 5 4 3 2 1 0]
895
sage: G
896
[ 0 1 2 3 4 5 6 -1]
897
[ 1 0 3 2 5 4 -1 -1]
898
[ 2 3 0 1 6 -1 4 -1]
899
[ 3 2 1 0 -1 -1 -1 -1]
900
[ 4 5 6 -1 0 1 2 -1]
901
[ 5 4 -1 -1 1 0 -1 -1]
902
[ 6 -1 4 -1 2 -1 0 -1]
903
[-1 -1 -1 -1 -1 -1 -1 -1]
904
"""
905
906
n = self.nrows()
907
908
from copy import copy
909
G = copy(self)
910
911
for r in range(n-1, -1, -1):
912
for c in range(n-1, -1, -1):
913
e = G[r, c]
914
G[r, c] = -1
915
916
if not G.dlxcpp_has_unique_completion():
917
G[r, c] = e
918
919
return G
920
921
def dlxcpp_has_unique_completion(self):
922
"""
923
Check if the partial latin square self of order n can be embedded
924
in precisely one latin square of order n.
925
926
EXAMPLES::
927
928
sage: from sage.combinat.matrices.latin import *
929
sage: back_circulant(2).dlxcpp_has_unique_completion()
930
True
931
sage: P = LatinSquare(2)
932
sage: P.dlxcpp_has_unique_completion()
933
False
934
sage: P[0, 0] = 0
935
sage: P.dlxcpp_has_unique_completion()
936
True
937
"""
938
return len(dlxcpp_find_completions(self, nr_to_find = 2)) == 1
939
940
def vals_in_row(self, r):
941
"""
942
Returns a dictionary with key e if and only if row r of self has
943
the symbol e.
944
945
EXAMPLES::
946
947
sage: from sage.combinat.matrices.latin import *
948
sage: B = back_circulant(3)
949
sage: B[0, 0] = -1
950
sage: back_circulant(3).vals_in_row(0)
951
{0: True, 1: True, 2: True}
952
"""
953
954
n = self.ncols()
955
vals_in_row = {}
956
957
for c in range(n):
958
e = self[r, c]
959
if e >= 0: vals_in_row[e] = True
960
961
return vals_in_row
962
963
def vals_in_col(self, c):
964
"""
965
Returns a dictionary with key e if and only if column c of self has
966
the symbol e.
967
968
EXAMPLES::
969
970
sage: from sage.combinat.matrices.latin import *
971
sage: B = back_circulant(3)
972
sage: B[0, 0] = -1
973
sage: back_circulant(3).vals_in_col(0)
974
{0: True, 1: True, 2: True}
975
"""
976
n = self.nrows()
977
vals_in_col = {}
978
979
for r in range(n):
980
e = self[r, c]
981
if e >= 0: vals_in_col[e] = True
982
983
return vals_in_col
984
985
def latex(self):
986
r"""
987
Returns LaTeX code for the latin square.
988
989
EXAMPLES::
990
991
sage: from sage.combinat.matrices.latin import *
992
sage: print back_circulant(3).latex()
993
\begin{array}{|c|c|c|}\hline 0 & 1 & 2\\\hline 1 & 2 & 0\\\hline 2 & 0 & 1\\\hline\end{array}
994
"""
995
996
a = ""
997
a += r"\begin{array}{" + self.ncols()*"|c" + "|}"
998
for r in range(self.nrows()):
999
a += r"\hline "
1000
for c in range(self.ncols()):
1001
s = self[r, c]
1002
1003
if s < 0: a += "~"
1004
else: a += str(s)
1005
1006
if c < self.ncols()-1: a += " & "
1007
else: a += "\\\\"
1008
a += r"\hline"
1009
a += r"\end{array}"
1010
return a
1011
1012
def disjoint_mate_dlxcpp_rows_and_map(self, allow_subtrade):
1013
"""
1014
Internal function for find_disjoint_mates.
1015
1016
EXAMPLES::
1017
1018
sage: from sage.combinat.matrices.latin import *
1019
sage: B = back_circulant(4)
1020
sage: B.disjoint_mate_dlxcpp_rows_and_map(allow_subtrade = True)
1021
([[0, 16, 32], [1, 17, 32], [2, 18, 32], [3, 19, 32], [4, 16, 33], [5, 17, 33], [6, 18, 33], [7, 19, 33], [8, 16, 34], [9, 17, 34], [10, 18, 34], [11, 19, 34], [12, 16, 35], [13, 17, 35], [14, 18, 35], [15, 19, 35], [0, 20, 36], [1, 21, 36], [2, 22, 36], [3, 23, 36], [4, 20, 37], [5, 21, 37], [6, 22, 37], [7, 23, 37], [8, 20, 38], [9, 21, 38], [10, 22, 38], [11, 23, 38], [12, 20, 39], [13, 21, 39], [14, 22, 39], [15, 23, 39], [0, 24, 40], [1, 25, 40], [2, 26, 40], [3, 27, 40], [4, 24, 41], [5, 25, 41], [6, 26, 41], [7, 27, 41], [8, 24, 42], [9, 25, 42], [10, 26, 42], [11, 27, 42], [12, 24, 43], [13, 25, 43], [14, 26, 43], [15, 27, 43], [0, 28, 44], [1, 29, 44], [2, 30, 44], [3, 31, 44], [4, 28, 45], [5, 29, 45], [6, 30, 45], [7, 31, 45], [8, 28, 46], [9, 29, 46], [10, 30, 46], [11, 31, 46], [12, 28, 47], [13, 29, 47], [14, 30, 47], [15, 31, 47]], {(9, 29, 46): (3, 2, 1), (13, 17, 35): (0, 3, 1), (7, 19, 33): (0, 1, 3), (14, 26, 43): (2, 3, 2), (0, 28, 44): (3, 0, 0), (5, 25, 41): (2, 1, 1), (11, 31, 46): (3, 2, 3), (14, 18, 35): (0, 3, 2), (11, 23, 38): (1, 2, 3), (5, 29, 45): (3, 1, 1), (13, 21, 39): (1, 3, 1), (1, 29, 44): (3, 0, 1), (0, 20, 36): (1, 0, 0), (12, 24, 43): (2, 3, 0), (8, 28, 46): (3, 2, 0), (12, 20, 39): (1, 3, 0), (11, 27, 42): (2, 2, 3), (6, 22, 37): (1, 1, 2), (1, 17, 32): (0, 0, 1), (10, 18, 34): (0, 2, 2), (12, 28, 47): (3, 3, 0), (1, 25, 40): (2, 0, 1), (10, 22, 38): (1, 2, 2), (5, 17, 33): (0, 1, 1), (3, 23, 36): (1, 0, 3), (6, 26, 41): (2, 1, 2), (9, 25, 42): (2, 2, 1), (7, 31, 45): (3, 1, 3), (15, 27, 43): (2, 3, 3), (3, 31, 44): (3, 0, 3), (8, 20, 38): (1, 2, 0), (2, 22, 36): (1, 0, 2), (3, 19, 32): (0, 0, 3), (9, 17, 34): (0, 2, 1), (15, 31, 47): (3, 3, 3), (8, 16, 34): (0, 2, 0), (14, 22, 39): (1, 3, 2), (4, 16, 33): (0, 1, 0), (14, 30, 47): (3, 3, 2), (2, 30, 44): (3, 0, 2), (4, 20, 37): (1, 1, 0), (6, 30, 45): (3, 1, 2), (12, 16, 35): (0, 3, 0), (15, 19, 35): (0, 3, 3), (5, 21, 37): (1, 1, 1), (4, 24, 41): (2, 1, 0), (13, 25, 43): (2, 3, 1), (0, 16, 32): (0, 0, 0), (15, 23, 39): (1, 3, 3), (7, 23, 37): (1, 1, 3), (6, 18, 33): (0, 1, 2), (10, 30, 46): (3, 2, 2), (13, 29, 47): (3, 3, 1), (11, 19, 34): (0, 2, 3), (1, 21, 36): (1, 0, 1), (7, 27, 41): (2, 1, 3), (0, 24, 40): (2, 0, 0), (10, 26, 42): (2, 2, 2), (3, 27, 40): (2, 0, 3), (2, 26, 40): (2, 0, 2), (9, 21, 38): (1, 2, 1), (8, 24, 42): (2, 2, 0), (4, 28, 45): (3, 1, 0), (2, 18, 32): (0, 0, 2)})
1022
"""
1023
1024
assert self.nrows() == self.ncols()
1025
1026
n = self.nrows()
1027
1028
# We will need 3n^2 columns in total:
1029
#
1030
# n^2 for the xCy columns
1031
# n^2 for the xRy columns
1032
# n^2 for the xy columns
1033
1034
dlx_rows = []
1035
cmap = {}
1036
1037
max_column_nr = -1
1038
1039
for r in range(n):
1040
valsrow = self.vals_in_row(r)
1041
1042
for c in range(n):
1043
valscol = self.vals_in_col(c)
1044
1045
# If this is an empty cell of self then we do nothing.
1046
if self[r, c] < 0: continue
1047
1048
for e in uniq(valsrow.keys() + valscol.keys()):
1049
# These should be constants
1050
c_OFFSET = e + c*n
1051
r_OFFSET = e + r*n + n*n
1052
xy_OFFSET = 2*n*n + r*n + c
1053
1054
cmap[(c_OFFSET, r_OFFSET, xy_OFFSET)] = (r,c,e)
1055
1056
# The disjoint mate has to be disjoint.
1057
if (not allow_subtrade) and self[r, c] == e: continue
1058
1059
# The permissible symbols must come from this row/column.
1060
if not(valsrow.has_key(e)): continue
1061
if not(valscol.has_key(e)): continue
1062
1063
dlx_rows.append([c_OFFSET, r_OFFSET, xy_OFFSET])
1064
1065
if max_column_nr < max(c_OFFSET, r_OFFSET, xy_OFFSET):
1066
max_column_nr = max(c_OFFSET, r_OFFSET, xy_OFFSET)
1067
1068
# We will have missed some columns. We
1069
# have to add 'dummy' rows so that the C++ DLX solver will find
1070
# a solution.
1071
used_columns = flatten(dlx_rows)
1072
for i in range(0, max_column_nr+1):
1073
if not i in used_columns:
1074
dlx_rows.append([i])
1075
1076
return dlx_rows, cmap
1077
1078
1079
def find_disjoint_mates(self, nr_to_find = None, allow_subtrade = False):
1080
r"""
1081
.. warning:::
1082
1083
If allow_subtrade is True then we may return a partial
1084
latin square that is *not* disjoint to self. In that case,
1085
use bitrade(P, Q) to get an actual bitrade.
1086
1087
EXAMPLES::
1088
1089
sage: from sage.combinat.matrices.latin import *
1090
sage: B = back_circulant(4)
1091
sage: g = B.find_disjoint_mates(allow_subtrade = True)
1092
sage: B1 = g.next()
1093
sage: B0, B1 = bitrade(B, B1)
1094
sage: assert is_bitrade(B0, B1)
1095
sage: print B0, "\n,\n", B1
1096
[-1 1 2 -1]
1097
[-1 2 -1 0]
1098
[-1 -1 -1 -1]
1099
[-1 0 1 2]
1100
,
1101
[-1 2 1 -1]
1102
[-1 0 -1 2]
1103
[-1 -1 -1 -1]
1104
[-1 1 2 0]
1105
"""
1106
1107
assert self.nrows() == self.ncols()
1108
1109
n = self.nrows()
1110
1111
dlx_rows, cmap = self.disjoint_mate_dlxcpp_rows_and_map(allow_subtrade)
1112
1113
nr_found = 0
1114
1115
for x in DLXCPP(dlx_rows):
1116
nr_found += 1
1117
1118
from copy import deepcopy
1119
Q = deepcopy(self)
1120
1121
for y in x:
1122
if len(dlx_rows[y]) == 1: continue # dummy row
1123
(r, c, e) = cmap[tuple(dlx_rows[y])]
1124
Q[r, c] = e
1125
1126
yield Q
1127
1128
if nr_to_find is not None and nr_found >= nr_to_find: return
1129
1130
def contained_in(self, Q):
1131
r"""
1132
Returns True if self is a subset of Q?
1133
1134
EXAMPLES::
1135
1136
sage: from sage.combinat.matrices.latin import *
1137
sage: P = elementary_abelian_2group(2)
1138
sage: P[0, 0] = -1
1139
sage: P.contained_in(elementary_abelian_2group(2))
1140
True
1141
sage: back_circulant(4).contained_in(elementary_abelian_2group(2))
1142
False
1143
"""
1144
1145
for r in range(self.nrows()):
1146
for c in range(self.ncols()):
1147
if self[r, c] >= 0 and Q[r, c] < 0: return False
1148
if self[r, c] >= 0 and (self[r, c] != Q[r, c]): return False
1149
return True
1150
1151
def genus(T1, T2):
1152
"""
1153
Returns the genus of hypermap embedding associated with the bitrade
1154
(T1, T2). Informally, we compute the [tau_1, tau_2, tau_3]
1155
permutation representation of the bitrade. Each cycle of tau_1,
1156
tau_2, and tau_3 gives a rotation scheme for a black, white, and
1157
star vertex (respectively). The genus then comes from Euler's
1158
formula. For more details see Carlo Hamalainen: Partitioning
1159
3-homogeneous latin bitrades. To appear in Geometriae Dedicata,
1160
available at http://arxiv.org/abs/0710.0938
1161
1162
EXAMPLES::
1163
1164
sage: from sage.combinat.matrices.latin import *
1165
sage: (a, b, c, G) = alternating_group_bitrade_generators(1)
1166
sage: (T1, T2) = bitrade_from_group(a, b, c, G)
1167
sage: genus(T1, T2)
1168
1
1169
sage: (a, b, c, G) = pq_group_bitrade_generators(3, 7)
1170
sage: (T1, T2) = bitrade_from_group(a, b, c, G)
1171
sage: genus(T1, T2)
1172
3
1173
"""
1174
1175
cells_map, t1, t2, t3 = tau123(T1, T2)
1176
return (len(t1.to_cycles()) + len(t2.to_cycles()) + len(t3.to_cycles()) - T1.nr_filled_cells() - 2)/(-2)
1177
1178
def tau123(T1, T2):
1179
"""
1180
Compute the tau_i representation for a bitrade (T1, T2). See the
1181
functions tau1, tau2, and tau3 for the mathematical definitions.
1182
1183
OUTPUT:
1184
1185
- (cells_map, t1, t2, t3)
1186
1187
where cells_map is a map to/from the filled cells of T1, and t1,
1188
t2, t3 are the tau1, tau2, tau3 permutations.
1189
1190
EXAMPLES::
1191
1192
sage: from sage.combinat.matrices.latin import *
1193
sage: (a, b, c, G) = pq_group_bitrade_generators(3, 7)
1194
sage: (T1, T2) = bitrade_from_group(a, b, c, G)
1195
sage: T1
1196
[ 0 1 2 -1 -1 -1 -1]
1197
[ 2 4 3 -1 -1 -1 -1]
1198
[ 4 0 5 -1 -1 -1 -1]
1199
[ 3 5 1 -1 -1 -1 -1]
1200
[ 6 3 0 -1 -1 -1 -1]
1201
[ 1 6 4 -1 -1 -1 -1]
1202
[ 5 2 6 -1 -1 -1 -1]
1203
sage: T2
1204
[ 2 0 1 -1 -1 -1 -1]
1205
[ 3 2 4 -1 -1 -1 -1]
1206
[ 5 4 0 -1 -1 -1 -1]
1207
[ 1 3 5 -1 -1 -1 -1]
1208
[ 0 6 3 -1 -1 -1 -1]
1209
[ 4 1 6 -1 -1 -1 -1]
1210
[ 6 5 2 -1 -1 -1 -1]
1211
sage: (cells_map, t1, t2, t3) = tau123(T1, T2)
1212
sage: cells_map
1213
{1: (0, 0), 2: (0, 1), 3: (0, 2), 4: (1, 0), 5: (1, 1), 6: (1, 2), 7: (2, 0), 8: (2, 1), 9: (2, 2), 10: (3, 0), 11: (3, 1), 12: (3, 2), 13: (4, 0), (2, 1): 8, 15: (4, 2), 16: (5, 0), 17: (5, 1), 18: (5, 2), 19: (6, 0), 20: (6, 1), 21: (6, 2), (5, 1): 17, (4, 0): 13, (1, 2): 6, (3, 0): 10, (5, 0): 16, (2, 2): 9, (4, 1): 14, (1, 1): 5, (3, 2): 12, (0, 0): 1, (6, 0): 19, 14: (4, 1), (4, 2): 15, (1, 0): 4, (0, 1): 2, (6, 1): 20, (3, 1): 11, (2, 0): 7, (6, 2): 21, (5, 2): 18, (0, 2): 3}
1214
sage: cells_map_as_square(cells_map, max(T1.nrows(), T1.ncols()))
1215
[ 1 2 3 -1 -1 -1 -1]
1216
[ 4 5 6 -1 -1 -1 -1]
1217
[ 7 8 9 -1 -1 -1 -1]
1218
[10 11 12 -1 -1 -1 -1]
1219
[13 14 15 -1 -1 -1 -1]
1220
[16 17 18 -1 -1 -1 -1]
1221
[19 20 21 -1 -1 -1 -1]
1222
sage: t1
1223
[2, 3, 1, 5, 6, 4, 8, 9, 7, 11, 12, 10, 14, 15, 13, 17, 18, 16, 20, 21, 19]
1224
sage: t2
1225
[4, 8, 12, 10, 20, 18, 19, 5, 15, 16, 14, 9, 1, 17, 6, 7, 2, 21, 13, 11, 3]
1226
sage: t3
1227
[15, 16, 20, 3, 7, 14, 18, 1, 11, 6, 19, 2, 21, 10, 8, 12, 13, 5, 9, 4, 17]
1228
1229
::
1230
1231
sage: t1.to_cycles()
1232
[(1, 2, 3), (4, 5, 6), (7, 8, 9), (10, 11, 12), (13, 14, 15), (16, 17, 18), (19, 20, 21)]
1233
sage: t2.to_cycles()
1234
[(1, 4, 10, 16, 7, 19, 13), (2, 8, 5, 20, 11, 14, 17), (3, 12, 9, 15, 6, 18, 21)]
1235
sage: t3.to_cycles()
1236
[(1, 15, 8), (2, 16, 12), (3, 20, 4), (5, 7, 18), (6, 14, 10), (9, 11, 19), (13, 21, 17)]
1237
1238
The product t1\*t2\*t3 is the identity, i.e. it fixes every point::
1239
1240
sage: len((t1*t2*t3).fixed_points()) == T1.nr_filled_cells()
1241
True
1242
"""
1243
1244
assert is_bitrade(T1, T2)
1245
1246
cells_map = T1.filled_cells_map()
1247
1248
t1 = tau1(T1, T2, cells_map)
1249
t2 = tau2(T1, T2, cells_map)
1250
t3 = tau3(T1, T2, cells_map)
1251
1252
return (cells_map, t1, t2, t3)
1253
1254
1255
def isotopism(p):
1256
"""
1257
Returns a Permutation object that represents an isotopism (for
1258
rows, columns or symbols of a partial latin square). Since matrices
1259
in Sage are indexed from 0, this function translates +1 to agree
1260
with the Permutation class. We also handle
1261
PermutationGroupElements.
1262
1263
EXAMPLES::
1264
1265
sage: from sage.combinat.matrices.latin import *
1266
sage: isotopism(5) # identity on 5 points
1267
[1, 2, 3, 4, 5]
1268
1269
::
1270
1271
sage: G = PermutationGroup(['(1,2,3)(4,5)'])
1272
sage: g = G.gen(0)
1273
sage: isotopism(g)
1274
[2, 3, 1, 5, 4]
1275
1276
::
1277
1278
sage: isotopism([0,3,2,1]) # 0 goes to 0, 1 goes to 3, etc.
1279
[1, 4, 3, 2]
1280
1281
::
1282
1283
sage: isotopism( (0,1,2) ) # single cycle, presented as a tuple
1284
[2, 3, 1]
1285
1286
::
1287
1288
sage: x = isotopism( ((0,1,2), (3,4)) ) # tuple of cycles
1289
sage: x
1290
[2, 3, 1, 5, 4]
1291
sage: x.to_cycles()
1292
[(1, 2, 3), (4, 5)]
1293
"""
1294
1295
# Identity isotopism on p points:
1296
if type(p) == Integer or type(p) == int:
1297
return Permutation(range(1, p+1))
1298
1299
if type(p) == PermutationGroupElement:
1300
# fixme Ask the Sage mailing list about the tuple/list issue!
1301
return Permutation(list(p.tuple()))
1302
1303
if type(p) == list:
1304
# We expect a list like [0,3,2,1] which means
1305
# that 0 goes to 0, 1 goes to 3, etc.
1306
return Permutation(map(lambda x: x+1, p))
1307
1308
if type(p) == tuple:
1309
# We have a single cycle:
1310
if type(p[0]) == Integer:
1311
return Permutation(tuple(map(lambda x: x+1, p)))
1312
1313
# We have a tuple of cycles:
1314
if type(p[0]) == tuple:
1315
x = isotopism(p[0])
1316
1317
for i in range(1, len(p)):
1318
x = x * isotopism(p[i])
1319
1320
return x
1321
1322
# Not sure what we got!
1323
raise NotImplemented
1324
1325
def cells_map_as_square(cells_map, n):
1326
"""
1327
Returns a LatinSquare with cells numbered from 1, 2, ... to given
1328
the dictionary cells_map.
1329
1330
.. note::
1331
1332
The value n should be the maximum of the number of rows and
1333
columns of the original partial latin square
1334
1335
EXAMPLES::
1336
1337
sage: from sage.combinat.matrices.latin import *
1338
sage: (a, b, c, G) = alternating_group_bitrade_generators(1)
1339
sage: (T1, T2) = bitrade_from_group(a, b, c, G)
1340
sage: T1
1341
[ 0 2 -1 1]
1342
[-1 0 1 3]
1343
[ 3 -1 0 2]
1344
[ 1 3 2 -1]
1345
1346
There are 12 filled cells in T::
1347
1348
sage: cells_map_as_square(T1.filled_cells_map(), max(T1.nrows(), T1.ncols()))
1349
[ 1 2 -1 3]
1350
[-1 4 5 6]
1351
[ 7 -1 8 9]
1352
[10 11 12 -1]
1353
"""
1354
1355
assert n > 1
1356
1357
L = LatinSquare(n, n)
1358
1359
for r in range(n):
1360
for c in range(n):
1361
try:
1362
L[r, c] = cells_map[ (r,c) ]
1363
except KeyError:
1364
# There is no cell (r,c) so skip it
1365
L[r, c] = -1
1366
1367
return L
1368
1369
def beta1(rce, T1, T2):
1370
"""
1371
Find the unique (x, c, e) in T2 such that (r, c, e) is in T1.
1372
1373
INPUT:
1374
1375
1376
- ``rce`` - tuple (or list) (r, c, e) in T1
1377
1378
- ``T1, T2`` - latin bitrade
1379
1380
1381
OUTPUT: (x, c, e) in T2.
1382
1383
EXAMPLES::
1384
1385
sage: from sage.combinat.matrices.latin import *
1386
sage: T1 = back_circulant(5)
1387
sage: x = isotopism( (0,1,2,3,4) )
1388
sage: y = isotopism(5) # identity
1389
sage: z = isotopism(5) # identity
1390
sage: T2 = T1.apply_isotopism(x, y, z)
1391
sage: is_bitrade(T1, T2)
1392
True
1393
sage: beta1([0, 0, 0], T1, T2)
1394
(1, 0, 0)
1395
"""
1396
1397
r = rce[0]
1398
c = rce[1]
1399
e = rce[2]
1400
1401
assert T1[r, c] == e
1402
assert e >= 0
1403
1404
for x in range(T1.nrows()):
1405
if T2[x, c] == e: return (x, c, e)
1406
1407
raise PairNotBitrade
1408
1409
def beta2(rce, T1, T2):
1410
"""
1411
Find the unique (r, x, e) in T2 such that (r, c, e) is in T1.
1412
1413
INPUT:
1414
1415
- ``rce`` - tuple (or list) (r, c, e) in T1
1416
1417
- ``T1, T2`` - latin bitrade
1418
1419
1420
OUTPUT:
1421
1422
- (r, x, e) in T2.
1423
1424
EXAMPLES::
1425
1426
sage: from sage.combinat.matrices.latin import *
1427
sage: T1 = back_circulant(5)
1428
sage: x = isotopism( (0,1,2,3,4) )
1429
sage: y = isotopism(5) # identity
1430
sage: z = isotopism(5) # identity
1431
sage: T2 = T1.apply_isotopism(x, y, z)
1432
sage: is_bitrade(T1, T2)
1433
True
1434
sage: beta2([0, 0, 0], T1, T2)
1435
(0, 1, 0)
1436
"""
1437
1438
r = rce[0]
1439
c = rce[1]
1440
e = rce[2]
1441
1442
assert T1[r, c] == e
1443
assert e >= 0
1444
1445
for x in range(T1.ncols()):
1446
if T2[r, x] == e: return (r, x, e)
1447
1448
raise PairNotBitrade
1449
1450
def beta3(rce, T1, T2):
1451
"""
1452
Find the unique (r, c, x) in T2 such that (r, c, e) is in T1.
1453
1454
INPUT:
1455
1456
1457
- ``rce`` - tuple (or list) (r, c, e) in T1
1458
1459
- ``T1, T2`` - latin bitrade
1460
1461
1462
OUTPUT:
1463
1464
- (r, c, x) in T2.
1465
1466
EXAMPLES::
1467
1468
sage: from sage.combinat.matrices.latin import *
1469
sage: T1 = back_circulant(5)
1470
sage: x = isotopism( (0,1,2,3,4) )
1471
sage: y = isotopism(5) # identity
1472
sage: z = isotopism(5) # identity
1473
sage: T2 = T1.apply_isotopism(x, y, z)
1474
sage: is_bitrade(T1, T2)
1475
True
1476
sage: beta3([0, 0, 0], T1, T2)
1477
(0, 0, 4)
1478
"""
1479
1480
r = rce[0]
1481
c = rce[1]
1482
e = rce[2]
1483
1484
assert T1[r, c] == e
1485
assert e >= 0
1486
1487
# fixme this range could be iffy if we
1488
# work with latin bitrade rectangles...
1489
for x in range(T1.nrows()):
1490
if T2[r, c] == x: return (r, c, x)
1491
1492
raise PairNotBitrade
1493
1494
def tau1(T1, T2, cells_map):
1495
r"""
1496
The definition of `\tau_1` is
1497
1498
.. math::
1499
1500
\tau_1 : T1 \rightarrow T1 \\
1501
\tau_1 = \beta_2^{-1} \beta_3
1502
1503
where the composition is left to right and `\beta_i : T2 \rightarrow T1`
1504
changes just the `i^{th}` coordinate of a triple.
1505
1506
EXAMPLES::
1507
1508
sage: from sage.combinat.matrices.latin import *
1509
sage: T1 = back_circulant(5)
1510
sage: x = isotopism( (0,1,2,3,4) )
1511
sage: y = isotopism(5) # identity
1512
sage: z = isotopism(5) # identity
1513
sage: T2 = T1.apply_isotopism(x, y, z)
1514
sage: is_bitrade(T1, T2)
1515
True
1516
sage: (cells_map, t1, t2, t3) = tau123(T1, T2)
1517
sage: t1 = tau1(T1, T2, cells_map)
1518
sage: t1
1519
[2, 3, 4, 5, 1, 7, 8, 9, 10, 6, 12, 13, 14, 15, 11, 17, 18, 19, 20, 16, 22, 23, 24, 25, 21]
1520
sage: t1.to_cycles()
1521
[(1, 2, 3, 4, 5), (6, 7, 8, 9, 10), (11, 12, 13, 14, 15), (16, 17, 18, 19, 20), (21, 22, 23, 24, 25)]
1522
"""
1523
1524
# The cells_map has both directions, i.e. integer to
1525
# cell and cell to integer, so the size of T1 is
1526
# just half of len(cells_map).
1527
x = (int(len(cells_map)/2) + 1) * [-1]
1528
1529
for r in range(T1.nrows()):
1530
for c in range(T1.ncols()):
1531
e = T1[r, c]
1532
1533
if e < 0: continue
1534
1535
(r2, c2, e2) = beta2( (r,c,e), T1, T2)
1536
(r3, c3, e3) = beta3( (r2,c2,e2), T2, T1)
1537
1538
x[ cells_map[(r,c)] ] = cells_map[ (r3,c3) ]
1539
1540
x.pop(0) # remove the head of the list since we
1541
# have permutations on 1..(something).
1542
1543
return Permutation(x)
1544
1545
def tau2(T1, T2, cells_map):
1546
r"""
1547
The definition of `\tau_2` is
1548
1549
.. math::
1550
1551
\tau_2 : T1 \rightarrow T1 \\
1552
\tau_2 = \beta_3^{-1} \beta_1
1553
1554
where the composition is left to right and `\beta_i : T2 \rightarrow T1`
1555
changes just the `i^{th}` coordinate of a triple.
1556
1557
EXAMPLES::
1558
1559
sage: from sage.combinat.matrices.latin import *
1560
sage: T1 = back_circulant(5)
1561
sage: x = isotopism( (0,1,2,3,4) )
1562
sage: y = isotopism(5) # identity
1563
sage: z = isotopism(5) # identity
1564
sage: T2 = T1.apply_isotopism(x, y, z)
1565
sage: is_bitrade(T1, T2)
1566
True
1567
sage: (cells_map, t1, t2, t3) = tau123(T1, T2)
1568
sage: t2 = tau2(T1, T2, cells_map)
1569
sage: t2
1570
[21, 22, 23, 24, 25, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]
1571
sage: t2.to_cycles()
1572
[(1, 21, 16, 11, 6), (2, 22, 17, 12, 7), (3, 23, 18, 13, 8), (4, 24, 19, 14, 9), (5, 25, 20, 15, 10)]
1573
"""
1574
1575
# The cells_map has both directions, i.e. integer to
1576
# cell and cell to integer, so the size of T1 is
1577
# just half of len(cells_map).
1578
x = (int(len(cells_map)/2) + 1) * [-1]
1579
1580
for r in range(T1.nrows()):
1581
for c in range(T1.ncols()):
1582
e = T1[r, c]
1583
1584
if e < 0: continue
1585
1586
(r2, c2, e2) = beta3( (r,c,e), T1, T2)
1587
(r3, c3, e3) = beta1( (r2,c2,e2), T2, T1)
1588
1589
x[ cells_map[(r,c)] ] = cells_map[ (r3,c3) ]
1590
1591
x.pop(0) # remove the head of the list since we
1592
# have permutations on 1..(something).
1593
1594
return Permutation(x)
1595
1596
def tau3(T1, T2, cells_map):
1597
r"""
1598
The definition of `\tau_3` is
1599
1600
.. math::
1601
1602
\tau_3 : T1 \rightarrow T1 \\
1603
\tau_3 = \beta_1^{-1} \beta_2
1604
1605
where the composition is left to right and `\beta_i : T2 \rightarrow T1`
1606
changes just the `i^{th}` coordinate of a triple.
1607
1608
EXAMPLES::
1609
1610
sage: from sage.combinat.matrices.latin import *
1611
sage: T1 = back_circulant(5)
1612
sage: x = isotopism( (0,1,2,3,4) )
1613
sage: y = isotopism(5) # identity
1614
sage: z = isotopism(5) # identity
1615
sage: T2 = T1.apply_isotopism(x, y, z)
1616
sage: is_bitrade(T1, T2)
1617
True
1618
sage: (cells_map, t1, t2, t3) = tau123(T1, T2)
1619
sage: t3 = tau3(T1, T2, cells_map)
1620
sage: t3
1621
[10, 6, 7, 8, 9, 15, 11, 12, 13, 14, 20, 16, 17, 18, 19, 25, 21, 22, 23, 24, 5, 1, 2, 3, 4]
1622
sage: t3.to_cycles()
1623
[(1, 10, 14, 18, 22), (2, 6, 15, 19, 23), (3, 7, 11, 20, 24), (4, 8, 12, 16, 25), (5, 9, 13, 17, 21)]
1624
"""
1625
1626
# The cells_map has both directions, i.e. integer to
1627
# cell and cell to integer, so the size of T1 is
1628
# just half of len(cells_map).
1629
x = (int(len(cells_map)/2) + 1) * [-1]
1630
1631
for r in range(T1.nrows()):
1632
for c in range(T1.ncols()):
1633
e = T1[r, c]
1634
1635
if e < 0: continue
1636
1637
(r2, c2, e2) = beta1( (r,c,e), T1, T2)
1638
(r3, c3, e3) = beta2( (r2,c2,e2), T2, T1)
1639
1640
x[ cells_map[(r,c)] ] = cells_map[ (r3,c3) ]
1641
1642
x.pop(0) # remove the head of the list since we
1643
# have permutations on 1..(something).
1644
1645
return Permutation(x)
1646
1647
def back_circulant(n):
1648
"""
1649
The back-circulant latin square of order n is the Cayley table for
1650
(Z_n, +), the integers under addition modulo n.
1651
1652
INPUT:
1653
1654
1655
- ``n`` - int; order of the latin square.
1656
1657
1658
EXAMPLES::
1659
1660
sage: from sage.combinat.matrices.latin import *
1661
sage: back_circulant(5)
1662
[0 1 2 3 4]
1663
[1 2 3 4 0]
1664
[2 3 4 0 1]
1665
[3 4 0 1 2]
1666
[4 0 1 2 3]
1667
"""
1668
1669
assert n >= 1
1670
1671
L = LatinSquare(n, n)
1672
1673
for r in range(n):
1674
for c in range(n):
1675
L[r, c] = (r + c) % n
1676
1677
return L
1678
1679
def forward_circulant(n):
1680
"""
1681
The forward-circulant latin square of order n is the Cayley table
1682
for the operation r + c = (n-c+r) mod n.
1683
1684
INPUT:
1685
1686
1687
- ``n`` - int; order of the latin square.
1688
1689
1690
EXAMPLES::
1691
1692
sage: from sage.combinat.matrices.latin import *
1693
sage: forward_circulant(5)
1694
[0 4 3 2 1]
1695
[1 0 4 3 2]
1696
[2 1 0 4 3]
1697
[3 2 1 0 4]
1698
[4 3 2 1 0]
1699
"""
1700
1701
assert n >= 1
1702
1703
L = LatinSquare(n, n)
1704
1705
for r in range(n):
1706
for c in range(n):
1707
L[r, c] = (n-c+r) % n
1708
1709
return L
1710
1711
def direct_product(L1, L2, L3, L4):
1712
"""
1713
The 'direct product' of four latin squares L1, L2, L3, L4 of order
1714
n is the latin square of order 2n consisting of
1715
1716
::
1717
1718
-----------
1719
| L1 | L2 |
1720
-----------
1721
| L3 | L4 |
1722
-----------
1723
1724
where the subsquares L2 and L3 have entries offset by n.
1725
1726
EXAMPLES::
1727
1728
sage: from sage.combinat.matrices.latin import *
1729
sage: direct_product(back_circulant(4), back_circulant(4), elementary_abelian_2group(2), elementary_abelian_2group(2))
1730
[0 1 2 3 4 5 6 7]
1731
[1 2 3 0 5 6 7 4]
1732
[2 3 0 1 6 7 4 5]
1733
[3 0 1 2 7 4 5 6]
1734
[4 5 6 7 0 1 2 3]
1735
[5 4 7 6 1 0 3 2]
1736
[6 7 4 5 2 3 0 1]
1737
[7 6 5 4 3 2 1 0]
1738
"""
1739
1740
assert L1.nrows() == L2.nrows() == L3.nrows() == L4.nrows()
1741
assert L1.ncols() == L2.ncols() == L3.ncols() == L4.ncols()
1742
assert L1.nrows() == L1.ncols()
1743
1744
n = L1.nrows()
1745
1746
D = LatinSquare(2*n, 2*n)
1747
1748
for r in range(n):
1749
for c in range(n):
1750
D[r, c] = L1[r, c]
1751
D[r, c+n] = L2[r, c] + n
1752
D[r+n, c] = L3[r, c] + n
1753
D[r+n, c+n] = L4[r, c]
1754
1755
return D
1756
1757
def elementary_abelian_2group(s):
1758
"""
1759
Returns the latin square based on the Cayley table for the
1760
elementary abelian 2-group of order 2s.
1761
1762
INPUT:
1763
1764
1765
- ``s`` - int; order of the latin square will be 2s.
1766
1767
1768
EXAMPLES::
1769
1770
sage: from sage.combinat.matrices.latin import *
1771
sage: elementary_abelian_2group(3)
1772
[0 1 2 3 4 5 6 7]
1773
[1 0 3 2 5 4 7 6]
1774
[2 3 0 1 6 7 4 5]
1775
[3 2 1 0 7 6 5 4]
1776
[4 5 6 7 0 1 2 3]
1777
[5 4 7 6 1 0 3 2]
1778
[6 7 4 5 2 3 0 1]
1779
[7 6 5 4 3 2 1 0]
1780
"""
1781
1782
assert s > 0
1783
1784
if s == 1:
1785
L = LatinSquare(2, 2)
1786
L[0, 0] = 0
1787
L[0, 1] = 1
1788
L[1, 0] = 1
1789
L[1, 1] = 0
1790
1791
return L
1792
else:
1793
L_prev = elementary_abelian_2group(s-1)
1794
L = LatinSquare(2**s, 2**s)
1795
1796
offset = L.nrows()/2
1797
1798
for r in range(L_prev.nrows()):
1799
for c in range(L_prev.ncols()):
1800
L[r, c] = L_prev[r, c]
1801
L[r+offset, c] = L_prev[r, c] + offset
1802
L[r, c+offset] = L_prev[r, c] + offset
1803
L[r+offset, c+offset] = L_prev[r, c]
1804
return L
1805
1806
def coin():
1807
"""
1808
Simulates a fair coin (returns True or False) using
1809
ZZ.random_element(2).
1810
1811
EXAMPLES::
1812
1813
sage: from sage.combinat.matrices.latin import *
1814
sage: x = coin()
1815
sage: x == 0 or x == 1
1816
True
1817
"""
1818
1819
return ZZ.random_element(2) == 0
1820
1821
1822
def next_conjugate(L):
1823
"""
1824
Permutes L[r, c] = e to the conjugate L[c, e] = r.
1825
1826
We assume that L is an n by n matrix and has values in the range 0,
1827
1, ..., n-1.
1828
1829
EXAMPLES::
1830
1831
sage: from sage.combinat.matrices.latin import *
1832
sage: L = back_circulant(6)
1833
sage: L
1834
[0 1 2 3 4 5]
1835
[1 2 3 4 5 0]
1836
[2 3 4 5 0 1]
1837
[3 4 5 0 1 2]
1838
[4 5 0 1 2 3]
1839
[5 0 1 2 3 4]
1840
sage: next_conjugate(L)
1841
[0 1 2 3 4 5]
1842
[5 0 1 2 3 4]
1843
[4 5 0 1 2 3]
1844
[3 4 5 0 1 2]
1845
[2 3 4 5 0 1]
1846
[1 2 3 4 5 0]
1847
sage: L == next_conjugate(next_conjugate(next_conjugate(L)))
1848
True
1849
"""
1850
1851
assert L.nrows() == L.ncols()
1852
1853
n = L.nrows()
1854
1855
C = LatinSquare(n, n)
1856
1857
for r in range(n):
1858
for c in range(n):
1859
e = L[r, c]
1860
assert e >= 0 and e < n
1861
1862
C[c, e] = r
1863
1864
return C
1865
1866
def row_containing_sym(L, c, x):
1867
"""
1868
Given an improper latin square L with L[r1, c] = L[r2, c] = x,
1869
return r1 or r2 with equal probability. This is an internal
1870
function and should only be used in LatinSquare_generator().
1871
1872
EXAMPLES::
1873
1874
sage: from sage.combinat.matrices.latin import *
1875
sage: L = matrix([(0, 1, 0, 3), (3, 0, 2, 1), (1, 0, 3, 2), (2, 3, 1, 0)])
1876
sage: L
1877
[0 1 0 3]
1878
[3 0 2 1]
1879
[1 0 3 2]
1880
[2 3 1 0]
1881
sage: c = row_containing_sym(L, 1, 0)
1882
sage: c == 1 or c == 2
1883
True
1884
"""
1885
1886
r1 = -1
1887
r2 = -1
1888
1889
for r in range(L.nrows()):
1890
if r1 >= 0 and r2 >= 0: break
1891
1892
if L[r, c] == x and r1 < 0:
1893
r1 = r
1894
continue
1895
1896
if L[r, c] == x and r2 < 0:
1897
r2 = r
1898
break
1899
1900
assert r1 >= 0 and r2 >= 0
1901
1902
if coin(): return r1
1903
else: return r2
1904
1905
def column_containing_sym(L, r, x):
1906
"""
1907
Given an improper latin square L with L[r, c1] = L[r, c2] = x,
1908
return c1 or c2 with equal probability. This is an internal
1909
function and should only be used in LatinSquare_generator().
1910
1911
EXAMPLES::
1912
1913
sage: from sage.combinat.matrices.latin import *
1914
sage: L = matrix([(1, 0, 2, 3), (0, 2, 3, 0), (2, 3, 0, 1), (3, 0, 1, 2)])
1915
sage: L
1916
[1 0 2 3]
1917
[0 2 3 0]
1918
[2 3 0 1]
1919
[3 0 1 2]
1920
sage: c = column_containing_sym(L, 1, 0)
1921
sage: c == 0 or c == 3
1922
True
1923
"""
1924
1925
c1 = -1
1926
c2 = -1
1927
1928
for c in range(L.ncols()):
1929
if c1 >= 0 and c2 >= 0: break
1930
1931
if L[r, c] == x and c1 < 0:
1932
c1 = c
1933
continue
1934
1935
if L[r, c] == x and c2 < 0:
1936
c2 = c
1937
break
1938
1939
assert c1 >= 0 and c2 >= 0
1940
1941
if coin(): return c1
1942
else: return c2
1943
1944
def LatinSquare_generator(L_start, check_assertions = False):
1945
"""
1946
Generator for a sequence of uniformly distributed latin squares,
1947
given L_start as the initial latin square. This code implements
1948
the Markov chain algorithm of Jacobson and Matthews (1996), see
1949
below for the BibTex entry. This generator will never throw the
1950
StopIteration exception, so it provides an infinite sequence of
1951
latin squares.
1952
1953
EXAMPLES:
1954
1955
Use the back circulant latin square of order 4 as the initial
1956
square and print the next two latin squares given by the Markov
1957
chain::
1958
1959
sage: from sage.combinat.matrices.latin import *
1960
sage: g = LatinSquare_generator(back_circulant(4))
1961
sage: g.next().is_latin_square()
1962
True
1963
1964
REFERENCE::
1965
1966
@article{MR1410617,
1967
AUTHOR = {Jacobson, Mark T. and Matthews, Peter},
1968
TITLE = {Generating uniformly distributed random {L}atin squares},
1969
JOURNAL = {J. Combin. Des.},
1970
FJOURNAL = {Journal of Combinatorial Designs},
1971
VOLUME = {4},
1972
YEAR = {1996},
1973
NUMBER = {6},
1974
PAGES = {405--437},
1975
ISSN = {1063-8539},
1976
MRCLASS = {05B15 (60J10)},
1977
MRNUMBER = {MR1410617 (98b:05021)},
1978
MRREVIEWER = {Lars D{\o}vling Andersen},
1979
}
1980
"""
1981
1982
if check_assertions: assert L_start.is_latin_square()
1983
1984
n = L_start.nrows()
1985
1986
r1 = r2 = c1 = c2 = x = y = z = -1
1987
proper = True
1988
1989
from copy import copy
1990
L = copy(L_start)
1991
1992
L_rce = L
1993
L_cer = LatinSquare(n, n)
1994
L_erc = LatinSquare(n, n)
1995
1996
while True:
1997
if proper:
1998
if check_assertions: assert L.is_latin_square()
1999
2000
#################################
2001
# Update the other two conjugates
2002
for r in range(n):
2003
for c in range(n):
2004
e = L[r, c]
2005
2006
L_cer[c, e] = r
2007
L_erc[e, r] = c
2008
#################################
2009
2010
yield L
2011
2012
r1 = ZZ.random_element(n)
2013
c1 = ZZ.random_element(n)
2014
x = L[r1, c1]
2015
2016
y = x
2017
while y == x:
2018
y = ZZ.random_element(n)
2019
2020
# Now find y in row r1 and column c1.
2021
if check_assertions:
2022
r2 = 0
2023
c2 = 0
2024
while L[r1, c2] != y: c2 += 1
2025
while L[r2, c1] != y: r2 += 1
2026
2027
assert L_erc[y, r1] == c2
2028
assert L_cer[c1, y] == r2
2029
2030
c2 = L_erc[y, r1]
2031
r2 = L_cer[c1, y]
2032
2033
if check_assertions: assert L[r1, c2] == y
2034
if check_assertions: assert L[r2, c1] == y
2035
2036
L[r1, c1] = y
2037
L[r1, c2] = x
2038
L[r2, c1] = x
2039
2040
# Now deal with the unknown point.
2041
# We want to form z + (y - x)
2042
z = L[r2, c2]
2043
2044
if z == x:
2045
L[r2, c2] = y
2046
else:
2047
# z and y have positive coefficients
2048
# x is the improper term with a negative coefficient
2049
proper = False
2050
else: # improper square,
2051
# L[r2, c2] = y + z - x
2052
# y and z are proper while x is the
2053
# improper symbol in the cell L[r2, c2].
2054
2055
r1 = row_containing_sym(L, c2, x)
2056
c1 = column_containing_sym(L, r2, x)
2057
2058
if check_assertions: assert L[r1, c2] == x
2059
if check_assertions: assert L[r2, c1] == x
2060
2061
# choose one of the proper symbols
2062
# uniformly at random (we will use whatever
2063
# lands in variable y).
2064
if coin():
2065
y, z = z, y
2066
2067
# Add/subtract the symbolic difference (y - x)
2068
L[r2, c2] = z
2069
L[r1, c2] = y
2070
L[r2, c1] = y
2071
2072
if L[r1, c1] == y:
2073
L[r1, c1] = x
2074
proper = True
2075
else: # got another improper square
2076
z = L[r1, c1]
2077
x, y = y, x
2078
r2 = r1
2079
c2 = c1
2080
2081
# Now we have L[r2, c2] = z+y-x as
2082
# usual
2083
proper = False # for emphasis
2084
2085
def group_to_LatinSquare(G):
2086
"""
2087
Construct a latin square on the symbols [0, 1, ..., n-1] for a
2088
group with an n by n Cayley table.
2089
2090
EXAMPLES::
2091
2092
sage: from sage.combinat.matrices.latin import *
2093
sage: group_to_LatinSquare(DihedralGroup(2))
2094
[0 1 2 3]
2095
[1 0 3 2]
2096
[2 3 0 1]
2097
[3 2 1 0]
2098
2099
::
2100
2101
sage: G = gap.Group(PermutationGroupElement((1,2,3)))
2102
sage: group_to_LatinSquare(G)
2103
[0 1 2]
2104
[1 2 0]
2105
[2 0 1]
2106
"""
2107
2108
if isinstance(G, GapElement):
2109
rows = map(lambda x: list(x), list(gap.MultiplicationTable(G)))
2110
new_rows = []
2111
2112
for x in rows:
2113
new_rows.append(map(lambda x: int(x)-1, x))
2114
2115
return matrix(new_rows)
2116
2117
# Otherwise we must have some kind of Sage permutation group object,
2118
# such as sage.groups.perm_gps.permgroup.PermutationGroup_generic
2119
# or maybe sage.groups.perm_gps.permgroup_named.
2120
2121
T = G.cayley_table()
2122
return matrix(ZZ, T.table())
2123
2124
def alternating_group_bitrade_generators(m):
2125
"""
2126
Construct generators a, b, c for the alternating group on 3m+1
2127
points, such that a\*b\*c = 1.
2128
2129
EXAMPLES::
2130
2131
sage: from sage.combinat.matrices.latin import *
2132
sage: a, b, c, G = alternating_group_bitrade_generators(1)
2133
sage: (a, b, c, G)
2134
((1,2,3), (1,4,2), (2,4,3), Permutation Group with generators [(1,2,3), (1,4,2)])
2135
sage: a*b*c
2136
()
2137
2138
::
2139
2140
sage: (T1, T2) = bitrade_from_group(a, b, c, G)
2141
sage: T1
2142
[ 0 2 -1 1]
2143
[-1 0 1 3]
2144
[ 3 -1 0 2]
2145
[ 1 3 2 -1]
2146
sage: T2
2147
[ 1 0 -1 2]
2148
[-1 3 0 1]
2149
[ 0 -1 2 3]
2150
[ 3 2 1 -1]
2151
"""
2152
2153
assert m >= 1
2154
2155
a = tuple(range(1, 2*m+1 + 1))
2156
2157
b = tuple(range(m+1, 0, -1) + range(2*m+2, 3*m+1 + 1))
2158
2159
a = PermutationGroupElement(a)
2160
b = PermutationGroupElement(b)
2161
c = PermutationGroupElement((a*b)**(-1))
2162
2163
G = PermutationGroup([a, b])
2164
2165
return (a, b, c, G)
2166
2167
def pq_group_bitrade_generators(p, q):
2168
"""
2169
Generators for a group of order pq where p and q are primes such
2170
that (q % p) == 1.
2171
2172
EXAMPLES::
2173
2174
sage: from sage.combinat.matrices.latin import *
2175
sage: pq_group_bitrade_generators(3,7)
2176
((2,3,5)(4,7,6), (1,2,3,4,5,6,7), (1,4,2)(3,5,6), Permutation Group with generators [(2,3,5)(4,7,6), (1,2,3,4,5,6,7)])
2177
"""
2178
2179
assert is_prime(p)
2180
assert is_prime(q)
2181
assert (q % p) == 1
2182
2183
# beta is a primitive root of the
2184
# congruence x^p = 1 mod q
2185
F = FiniteField(q)
2186
fgen = F.multiplicative_generator()
2187
beta = fgen**((q-1)/p)
2188
2189
assert beta != 1
2190
assert (beta**p % q) == 1
2191
2192
Q = tuple(range(1, q+1))
2193
2194
P = []
2195
seenValues = {}
2196
for i in range(2, q):
2197
if seenValues.has_key(i):
2198
continue
2199
2200
cycle = []
2201
for k in range(p):
2202
x = (1 + (i-1)*beta**k) % q
2203
if x == 0: x = q
2204
2205
seenValues[x] = True
2206
cycle.append(x)
2207
P.append(tuple(map(Integer,cycle)))
2208
2209
G = PermutationGroup([P, Q])
2210
assert G.order() == p*q
2211
assert not G.is_abelian()
2212
2213
a = PermutationGroupElement(P)
2214
b = PermutationGroupElement(Q)
2215
c = PermutationGroupElement((a*b)**(-1))
2216
2217
return (a, b, c, PermutationGroup([P, Q]))
2218
2219
def p3_group_bitrade_generators(p):
2220
"""
2221
Generators for a group of order p3 where p is a prime.
2222
2223
EXAMPLES::
2224
2225
sage: from sage.combinat.matrices.latin import *
2226
sage: p3_group_bitrade_generators(3)
2227
((2,6,7)(3,8,9), (1,2,3)(4,7,8)(5,6,9), (1,9,2)(3,7,4)(5,8,6), Permutation Group with generators [(2,6,7)(3,8,9), (1,2,3)(4,7,8)(5,6,9)])
2228
"""
2229
2230
assert is_prime(p)
2231
2232
F = gap.new("FreeGroup(3)")
2233
2234
a = F.gen(1)
2235
b = F.gen(2)
2236
c = F.gen(3)
2237
2238
rels = []
2239
rels.append( a**p )
2240
rels.append( b**p )
2241
rels.append( c**p )
2242
rels.append( a*b*((b*a*c)**(-1)) )
2243
rels.append( c*a*((a*c)**(-1)) )
2244
rels.append( c*b*((b*c)**(-1)) )
2245
2246
G = F.FactorGroupFpGroupByRels(rels)
2247
2248
iso = gap.IsomorphismPermGroup(G)
2249
2250
x = PermutationGroupElement(gap.Image(iso, G.gen(1)))
2251
y = PermutationGroupElement(gap.Image(iso, G.gen(2)))
2252
2253
return (x, y, (x*y)**(-1), PermutationGroup([x, y]))
2254
2255
def check_bitrade_generators(a, b, c):
2256
"""
2257
Three group elements a, b, c will generate a bitrade if a\*b\*c = 1
2258
and the subgroups a, b, c intersect (pairwise) in just the
2259
identity.
2260
2261
EXAMPLES::
2262
2263
sage: from sage.combinat.matrices.latin import *
2264
sage: a, b, c, G = p3_group_bitrade_generators(3)
2265
sage: check_bitrade_generators(a, b, c)
2266
True
2267
sage: check_bitrade_generators(a, b, gap('()'))
2268
False
2269
"""
2270
2271
A = PermutationGroup([a])
2272
B = PermutationGroup([b])
2273
C = PermutationGroup([c])
2274
2275
if a*b != c**(-1): return False
2276
2277
X = gap.Intersection(gap.Intersection(A, B), C)
2278
return X.Size() == 1
2279
2280
def is_bitrade(T1, T2):
2281
"""
2282
Combinatorially, a pair (T1, T2) of partial latin squares is a
2283
bitrade if they are disjoint, have the same shape, and have row and
2284
column balance. For definitions of each of these terms see the
2285
relevant function in this file.
2286
2287
EXAMPLES::
2288
2289
sage: from sage.combinat.matrices.latin import *
2290
sage: T1 = back_circulant(5)
2291
sage: x = isotopism( (0,1,2,3,4) )
2292
sage: y = isotopism(5) # identity
2293
sage: z = isotopism(5) # identity
2294
sage: T2 = T1.apply_isotopism(x, y, z)
2295
sage: is_bitrade(T1, T2)
2296
True
2297
"""
2298
2299
if not is_disjoint(T1, T2): return False
2300
if not is_same_shape(T1, T2): return False
2301
if not is_row_and_col_balanced(T1, T2): return False
2302
2303
return True
2304
2305
def is_primary_bitrade(a, b, c, G):
2306
"""
2307
A bitrade generated from elements a, b, c is primary if a, b, c =
2308
G.
2309
2310
EXAMPLES::
2311
2312
sage: from sage.combinat.matrices.latin import *
2313
sage: (a, b, c, G) = p3_group_bitrade_generators(5)
2314
sage: is_primary_bitrade(a, b, c, G)
2315
True
2316
"""
2317
2318
H = PermutationGroup([a, b, c])
2319
2320
return G == H
2321
2322
def tau_to_bitrade(t1, t2, t3):
2323
"""
2324
Given permutations t1, t2, t3 that represent a latin bitrade,
2325
convert them to an explicit latin bitrade (T1, T2). The result is
2326
unique up to isotopism.
2327
2328
EXAMPLE::
2329
2330
sage: from sage.combinat.matrices.latin import *
2331
sage: T1 = back_circulant(5)
2332
sage: x = isotopism( (0,1,2,3,4) )
2333
sage: y = isotopism(5) # identity
2334
sage: z = isotopism(5) # identity
2335
sage: T2 = T1.apply_isotopism(x, y, z)
2336
sage: _, t1, t2, t3 = tau123(T1, T2)
2337
sage: U1, U2 = tau_to_bitrade(t1, t2, t3)
2338
sage: assert is_bitrade(U1, U2)
2339
sage: U1
2340
[0 1 2 3 4]
2341
[1 2 3 4 0]
2342
[2 3 4 0 1]
2343
[3 4 0 1 2]
2344
[4 0 1 2 3]
2345
sage: U2
2346
[4 0 1 2 3]
2347
[0 1 2 3 4]
2348
[1 2 3 4 0]
2349
[2 3 4 0 1]
2350
[3 4 0 1 2]
2351
"""
2352
2353
c1 = t1.to_cycles()
2354
c2 = t2.to_cycles()
2355
c3 = t3.to_cycles()
2356
2357
pt_to_cycle1 = {}
2358
pt_to_cycle2 = {}
2359
pt_to_cycle3 = {}
2360
2361
for i in range(len(c1)):
2362
for j in range(len(c1[i])):
2363
pt_to_cycle1[c1[i][j]] = i
2364
2365
for i in range(len(c2)):
2366
for j in range(len(c2[i])):
2367
pt_to_cycle2[c2[i][j]] = i
2368
2369
for i in range(len(c3)):
2370
for j in range(len(c3[i])):
2371
pt_to_cycle3[c3[i][j]] = i
2372
2373
n = max(len(c1), len(c2), len(c3))
2374
2375
T1 = LatinSquare(n)
2376
T2 = LatinSquare(n)
2377
2378
for r in range(len(c1)):
2379
for c in range(len(c2)):
2380
for s in range(len(c3)):
2381
nr_common = len(reduce(set.intersection, \
2382
[set(c1[r]), set(c2[c]), set(c3[s])]))
2383
assert nr_common in [0, 1]
2384
2385
if nr_common == 1: T1[r, c] = s
2386
2387
for cycle in c1:
2388
for pt1 in cycle:
2389
pt2 = t1[pt1 - 1]
2390
pt3 = t2[pt2 - 1]
2391
assert t3[pt3 - 1] == pt1
2392
2393
r = pt_to_cycle1[pt1]
2394
c = pt_to_cycle2[pt2]
2395
s = pt_to_cycle3[pt3]
2396
2397
T2[r, c] = s
2398
2399
return T1, T2
2400
2401
2402
def bitrade_from_group(a, b, c, G):
2403
"""
2404
Given group elements a, b, c in G such that abc = 1 and the
2405
subgroups a, b, c intersect (pairwise) only in the identity,
2406
construct a bitrade (T1, T2) where rows, columns, and symbols
2407
correspond to cosets of a, b, and c, respectively.
2408
2409
EXAMPLES::
2410
2411
sage: from sage.combinat.matrices.latin import *
2412
sage: a, b, c, G = alternating_group_bitrade_generators(1)
2413
sage: (T1, T2) = bitrade_from_group(a, b, c, G)
2414
sage: T1
2415
[ 0 2 -1 1]
2416
[-1 0 1 3]
2417
[ 3 -1 0 2]
2418
[ 1 3 2 -1]
2419
sage: T2
2420
[ 1 0 -1 2]
2421
[-1 3 0 1]
2422
[ 0 -1 2 3]
2423
[ 3 2 1 -1]
2424
"""
2425
2426
hom = gap.ActionHomomorphism(G, gap.RightCosets(G, gap.TrivialSubgroup(G)), gap.OnRight)
2427
2428
t1 = gap.Image(hom, a)
2429
t2 = gap.Image(hom, b)
2430
t3 = gap.Image(hom, c)
2431
2432
t1 = Permutation(str(t1).replace('\n', ''))
2433
t2 = Permutation(str(t2).replace('\n', ''))
2434
t3 = Permutation(str(t3).replace('\n', ''))
2435
2436
return tau_to_bitrade(t1, t2, t3)
2437
2438
def is_disjoint(T1, T2):
2439
"""
2440
The partial latin squares T1 and T2 are disjoint if T1[r, c] !=
2441
T2[r, c] or T1[r, c] == T2[r, c] == -1 for each cell [r, c].
2442
2443
EXAMPLES::
2444
2445
sage: from sage.combinat.matrices.latin import is_disjoint, back_circulant, isotopism
2446
sage: is_disjoint(back_circulant(2), back_circulant(2))
2447
False
2448
2449
::
2450
2451
sage: T1 = back_circulant(5)
2452
sage: x = isotopism( (0,1,2,3,4) )
2453
sage: y = isotopism(5) # identity
2454
sage: z = isotopism(5) # identity
2455
sage: T2 = T1.apply_isotopism(x, y, z)
2456
sage: is_disjoint(T1, T2)
2457
True
2458
"""
2459
2460
for i in range(T1.nrows()):
2461
for j in range(T1.ncols()):
2462
if T1[i, j] < 0 and T2[i, j] < 0: continue
2463
2464
if T1[i, j] == T2[i, j]:
2465
return False
2466
2467
return True
2468
2469
def is_same_shape(T1, T2):
2470
"""
2471
Two partial latin squares T1, T2 have the same shape if T1[r, c] =
2472
0 if and only if T2[r, c] = 0.
2473
2474
EXAMPLES::
2475
2476
sage: from sage.combinat.matrices.latin import *
2477
sage: is_same_shape(elementary_abelian_2group(2), back_circulant(4))
2478
True
2479
sage: is_same_shape(LatinSquare(5), LatinSquare(5))
2480
True
2481
sage: is_same_shape(forward_circulant(5), LatinSquare(5))
2482
False
2483
"""
2484
2485
for i in range(T1.nrows()):
2486
for j in range(T1.ncols()):
2487
if T1[i, j] < 0 and T2[i, j] < 0: continue
2488
if T1[i, j] >= 0 and T2[i, j] >= 0: continue
2489
2490
return False
2491
2492
return True
2493
2494
def is_row_and_col_balanced(T1, T2):
2495
"""
2496
Partial latin squares T1 and T2 are balanced if the symbols
2497
appearing in row r of T1 are the same as the symbols appearing in
2498
row r of T2, for each r, and if the same condition holds on
2499
columns.
2500
2501
EXAMPLES::
2502
2503
sage: from sage.combinat.matrices.latin import *
2504
sage: T1 = matrix([[0,1,-1,-1], [-1,-1,-1,-1], [-1,-1,-1,-1], [-1,-1,-1,-1]])
2505
sage: T2 = matrix([[0,1,-1,-1], [-1,-1,-1,-1], [-1,-1,-1,-1], [-1,-1,-1,-1]])
2506
sage: is_row_and_col_balanced(T1, T2)
2507
True
2508
sage: T2 = matrix([[0,3,-1,-1], [-1,-1,-1,-1], [-1,-1,-1,-1], [-1,-1,-1,-1]])
2509
sage: is_row_and_col_balanced(T1, T2)
2510
False
2511
"""
2512
2513
for r in range(T1.nrows()):
2514
val1 = set(filter(lambda x: x >= 0, T1.row(r)))
2515
val2 = set(filter(lambda x: x >= 0, T2.row(r)))
2516
2517
if val1 != val2: return False
2518
2519
for c in range(T1.ncols()):
2520
val1 = set(filter(lambda x: x >= 0, T1.column(c)))
2521
val2 = set(filter(lambda x: x >= 0, T2.column(c)))
2522
2523
if val1 != val2: return False
2524
2525
return True
2526
2527
def dlxcpp_rows_and_map(P):
2528
"""
2529
Internal function for dlxcpp_find_completions. Given a partial
2530
latin square P we construct a list of rows of a 0-1 matrix M such
2531
that an exact cover of M corresponds to a completion of P to a
2532
latin square.
2533
2534
EXAMPLES::
2535
2536
sage: from sage.combinat.matrices.latin import *
2537
sage: dlxcpp_rows_and_map(LatinSquare(2))
2538
([[0, 4, 8], [1, 5, 8], [2, 4, 9], [3, 5, 9], [0, 6, 10], [1, 7, 10], [2, 6, 11], [3, 7, 11]], {(2, 4, 9): (0, 1, 0), (1, 5, 8): (0, 0, 1), (1, 7, 10): (1, 0, 1), (0, 6, 10): (1, 0, 0), (3, 7, 11): (1, 1, 1), (2, 6, 11): (1, 1, 0), (0, 4, 8): (0, 0, 0), (3, 5, 9): (0, 1, 1)})
2539
"""
2540
2541
assert P.nrows() == P.ncols()
2542
2543
n = P.nrows()
2544
2545
# We will need 3n^2 columns in total:
2546
#
2547
# n^2 for the xCy columns
2548
# n^2 for the xRy columns
2549
# n^2 for the xy columns
2550
2551
dlx_rows = []
2552
cmap = {}
2553
2554
for r in range(n):
2555
valsrow = P.vals_in_row(r)
2556
2557
for c in range(n):
2558
valscol = P.vals_in_col(c)
2559
2560
for e in range(n):
2561
# These should be constants
2562
c_OFFSET = e + c*n
2563
r_OFFSET = e + r*n + n*n
2564
xy_OFFSET = 2*n*n + r*n + c
2565
2566
cmap[(c_OFFSET, r_OFFSET, xy_OFFSET)] = (r,c,e)
2567
2568
#print "possibility: ", r, c, e, "offsets:", c_OFFSET, r_OFFSET, xy_OFFSET
2569
2570
#if P[r, c] >= 0: continue
2571
2572
# We only want the correct value to pop in here
2573
if P[r, c] >= 0 and P[r, c] != e: continue
2574
2575
if P[r, c] < 0 and valsrow.has_key(e): continue
2576
if P[r, c] < 0 and valscol.has_key(e): continue
2577
2578
dlx_rows.append([c_OFFSET, r_OFFSET, xy_OFFSET])
2579
2580
return dlx_rows, cmap
2581
2582
2583
def dlxcpp_find_completions(P, nr_to_find = None):
2584
"""
2585
Returns a list of all latin squares L of the same order as P such
2586
that P is contained in L. The optional parameter nr_to_find
2587
limits the number of latin squares that are found.
2588
2589
EXAMPLES::
2590
2591
sage: from sage.combinat.matrices.latin import *
2592
sage: dlxcpp_find_completions(LatinSquare(2))
2593
[[0 1]
2594
[1 0], [1 0]
2595
[0 1]]
2596
2597
::
2598
2599
sage: dlxcpp_find_completions(LatinSquare(2), 1)
2600
[[0 1]
2601
[1 0]]
2602
"""
2603
2604
2605
assert P.nrows() == P.ncols()
2606
2607
n = P.nrows()
2608
2609
dlx_rows, cmap = dlxcpp_rows_and_map(P)
2610
2611
SOLUTIONS = {}
2612
for x in DLXCPP(dlx_rows):
2613
x.sort()
2614
SOLUTIONS[tuple(x)] = True
2615
2616
if nr_to_find is not None and len(SOLUTIONS) >= nr_to_find: break
2617
2618
comps = []
2619
2620
for i in SOLUTIONS.keys():
2621
soln = list(i)
2622
2623
from copy import deepcopy
2624
Q = deepcopy(P)
2625
2626
for x in soln:
2627
(r, c, e) = cmap[tuple(dlx_rows[x])]
2628
2629
if Q[r, c] >= 0:
2630
assert Q[r, c] == e
2631
else:
2632
Q[r, c] = e
2633
2634
comps.append(Q)
2635
2636
return comps
2637
2638
def bitrade(T1, T2):
2639
r"""
2640
Form the bitrade (Q1, Q2) from (T1, T2) by setting empty the cells
2641
(r, c) such that T1[r, c] == T2[r, c].
2642
2643
EXAMPLES::
2644
2645
sage: from sage.combinat.matrices.latin import *
2646
sage: B1 = back_circulant(5)
2647
sage: alpha = isotopism((0,1,2,3,4))
2648
sage: beta = isotopism((1,0,2,3,4))
2649
sage: gamma = isotopism((2,1,0,3,4))
2650
sage: B2 = B1.apply_isotopism(alpha, beta, gamma)
2651
sage: T1, T2 = bitrade(B1, B2)
2652
sage: T1
2653
[ 0 1 -1 3 4]
2654
[ 1 -1 -1 4 0]
2655
[ 2 -1 4 0 1]
2656
[ 3 4 0 1 2]
2657
[ 4 0 1 2 3]
2658
sage: T2
2659
[ 3 4 -1 0 1]
2660
[ 0 -1 -1 1 4]
2661
[ 1 -1 0 4 2]
2662
[ 4 0 1 2 3]
2663
[ 2 1 4 3 0]
2664
"""
2665
assert T1.nrows() == T1.ncols()
2666
assert T2.nrows() == T2.ncols()
2667
assert T1.nrows() == T2.nrows()
2668
2669
n = T1.nrows()
2670
2671
from copy import copy
2672
Q1 = copy(T1)
2673
Q2 = copy(T2)
2674
2675
for r in range(n):
2676
for c in range(n):
2677
if T1[r, c] == T2[r, c]:
2678
Q1[r, c] = -1
2679
Q2[r, c] = -1
2680
2681
return Q1, Q2
2682
2683
2684
2685