trac #16231: Equivalence between OA/TD/MOLS

src/sage/combinat/designs/latin_squares.py

@@ -57,7 +57,7 @@ def are_mutually_orthogonal_latin_squares(l, verbose=False):
         matrices 0 and 2 are not orthogonal
         False
 
-        sage: m = designs.mutually_orthogonal_latin_squares(8)
+        sage: m = designs.mutually_orthogonal_latin_squares(8,7)
         sage: are_mutually_orthogonal_latin_squares(m)
         True
     """
@@ -91,8 +91,7 @@ def are_mutually_orthogonal_latin_squares(l, verbose=False):
 
     return True
 
-
-def mutually_orthogonal_latin_squares(n,k=None, partitions = False, availability=False):
+def mutually_orthogonal_latin_squares(n,k, partitions = False, check = True, availability=False, who_asked=tuple()):
     r"""
     Returns `k` Mutually Orthogonal `n\times n` Latin Squares (MOLS).
 
@@ -105,13 +104,7 @@ def mutually_orthogonal_latin_squares(n,k=None, partitions = False, availability
 
     - ``n`` (integer) -- size of the latin square.
 
-    - ``k`` (integer) -- returns `k` MOLS. If set to ``None`` (default), returns
-      the maximum number of MOLS that Sage can build.
-
-      .. WARNING::
-
-          This has no reason to be the maximum number of `n\times n` MOLS, just
-          the best Sage can do !
+    - ``k`` (integer) -- number of MOLS.
 
     - ``partition`` (boolean) -- a Latin Square can be seen as 3 partitions of
       the `n^2` cells of the array into `n` sets of size `n`, respectively :
@@ -133,9 +126,19 @@ def mutually_orthogonal_latin_squares(n,k=None, partitions = False, availability
       to build such a collection. This should be much faster than actually
       building it.
 
+    - ``check`` -- (boolean) Whether to check that output is correct before
+      returning it. As this is expected to be useless (but we are cautious
+      guys), you may want to disable it whenever you want speed. Set to
+      ``True`` by default.
+
+    - ``who_asked`` (internal use only) -- because of the equivalence between
+      OA/TD/MOLS, each of the three constructors calls the others. We must keep
+      track of who calls who in order to avoid infinite loops. ``who_asked`` is
+      the tuple of the other functions that were called before this one.
+
     EXAMPLES::
 
-        sage: designs.mutually_orthogonal_latin_squares(5)
+        sage: designs.mutually_orthogonal_latin_squares(5,4)
         [
         [0 1 2 3 4]  [0 1 2 3 4]  [0 1 2 3 4]  [0 1 2 3 4]
         [3 0 1 4 2]  [4 3 0 2 1]  [1 2 4 0 3]  [2 4 3 1 0]
@@ -184,22 +187,21 @@ def mutually_orthogonal_latin_squares(n,k=None, partitions = False, availability
     """
     from sage.rings.finite_rings.constructor import FiniteField
     from sage.combinat.designs.block_design import AffineGeometryDesign
+    from sage.combinat.designs.orthogonal_arrays import orthogonal_array
     from sage.rings.arith import is_prime_power
     from sage.matrix.constructor import Matrix
     from sage.rings.arith import factor
 
-    if k is not None and k >= n:
+    if k >= n:
         if availability:
             return False
         else:
             raise ValueError("There exist at most n-1 MOLS of size n.")
 
     if is_prime_power(n):
         if availability:
-            return n-1 if k is None else True
+            return True
 
-        if k is None:
-            k = n-1
         # Section 6.4.1 of [Stinson2004]
         Fp = FiniteField(n,'x')
         B = AffineGeometryDesign(2,1,Fp).blocks()
@@ -210,9 +212,6 @@ def mutually_orthogonal_latin_squares(n,k=None, partitions = False, availability
                     p.append(b)
                     break
 
-        if partitions:
-            return parallel_classes
-
         coord = {v:i
                  for i,L in enumerate(parallel_classes[0]) for v in L}
         coord = {v:(coord[v],i)
@@ -221,26 +220,26 @@ def mutually_orthogonal_latin_squares(n,k=None, partitions = False, availability
         matrices = []
         for P in parallel_classes[2:]:
             matrices.append(Matrix({coord[v]:i for i,L in enumerate(P) for v in L }))
-        return matrices
-    else:
-        # Theorem 6.33 of [Stinson2004], MacNeish's theorem.
-        subcases = [p**i for p,i in factor(n)]
-        s = min(subcases)-1
-        if k is None:
-            k = s
-            if availability:
-                return k
-        elif k > s:
-            if availability:
-                return False
-            else:
-                raise NotImplementedError("I don't know how to build these MOLS.")
-        elif availability:
+
+        if partitions:
+            partitions = parallel_classes
+
+    elif (orthogonal_array not in who_asked and
+        orthogonal_array(k+2,n,availability=True,who_asked = who_asked+(mutually_orthogonal_latin_squares,))):
+        if availability:
             return True
+        OA = orthogonal_array(k+2,n,check=False, who_asked = who_asked+(mutually_orthogonal_latin_squares,))
+        OA.sort() # make sure that the first two columns are "11, 12, ..., 1n, 21, 22, ..."
 
-        subcalls = [mutually_orthogonal_latin_squares(p,k) for p in subcases]
-        matrices = [latin_square_product(*[sc[i] for sc in subcalls])
-                    for i in range(k)]
+        # We first define matrices as lists of n^2 values
+        matrices = [[] for _ in range(k)]
+        for L in OA:
+            for i in range(2,k+2):
+                matrices[i-2].append(L[i])
+
+        # The real matrices
+        matrices = [[M[i*n:(i+1)*n] for i in range(n)] for M in matrices]
+        matrices = [Matrix(M) for M in matrices]
 
         if partitions:
             partitions = [[[i*n+j for j in range(n)] for i in range(n)],
@@ -252,10 +251,19 @@ def mutually_orthogonal_latin_squares(n,k=None, partitions = False, availability
                         partition[m[i,j]].append(i*n+j)
                 partitions.append(partition)
 
-            return partitions
-
+    else:
+        if availability:
+            return False
         else:
-            return matrices
+            raise NotImplementedError("I don't know how to build these MOLS!")
+
+    if check:
+        assert are_mutually_orthogonal_latin_squares(matrices)
+
+    if partitions:
+        return partitions
+    else:
+        return matrices
 
 def latin_square_product(M,N,*others):
     r"""
@@ -277,7 +285,7 @@ def latin_square_product(M,N,*others):
     EXAMPLES::
 
         sage: from sage.combinat.designs.latin_squares import latin_square_product
-        sage: m=designs.mutually_orthogonal_latin_squares(4)[0]
+        sage: m=designs.mutually_orthogonal_latin_squares(4,3)[0]
         sage: latin_square_product(m,m,m)
         64 x 64 sparse matrix over Integer Ring
     """

src/sage/combinat/designs/orthogonal_arrays.py

@@ -9,7 +9,7 @@
 """
 from sage.misc.cachefunc import cached_function
 
-def transversal_design(k,n,check=True,availability=False):
+def transversal_design(k,n,check=True,availability=False, who_asked=tuple()):
     r"""
     Return a transversal design of parameters `k,n`.
 
@@ -43,6 +43,11 @@ def transversal_design(k,n,check=True,availability=False):
       to build such a design. This should be much faster than actually building
       it.
 
+    - ``who_asked`` (internal use only) -- because of the equivalence between
+      OA/TD/MOLS, each of the three constructors calls the others. We must keep
+      track of who calls who in order to avoid infinite loops. ``who_asked`` is
+      the tuple of the other functions that were called before this one.
+
     .. NOTE::
 
         This function returns transversal designs with
@@ -75,18 +80,20 @@ def transversal_design(k,n,check=True,availability=False):
     """
     if n == 12 and k <= 6:
         TD = [l[:k] for l in TD6_12()]
-    # Section 6.6 of [Stinson2004]
-    elif orthogonal_array(k,n, check = False, availability = True):
-        if availability:
-            return True
-        OA = orthogonal_array(k,n, check = False)
-        TD = [[i*n+c for i,c in enumerate(l)] for l in OA]
 
     elif find_wilson_decomposition(k,n):
         if availability:
             return True
         TD = wilson_construction(*find_wilson_decomposition(k,n), check = False)
 
+    # Section 6.6 of [Stinson2004]
+    elif (orthogonal_array not in who_asked and
+          orthogonal_array(k, n, availability=True, who_asked = who_asked + (transversal_design,))):
+        if availability:
+            return True
+        OA = orthogonal_array(k,n, check = False, who_asked = who_asked + (transversal_design,))
+        TD = [[i*n+c for i,c in enumerate(l)] for l in OA]
+
     else:
         if availability:
             return False
@@ -327,7 +334,7 @@ def wilson_construction(k,m,t,u, check = True):
 
     return TD
 
-def orthogonal_array(k,n,t=2,check=True,availability=False):
+def orthogonal_array(k,n,t=2,check=True,availability=False,who_asked=tuple()):
     r"""
     Return an orthogonal array of parameters `k,n,t`.
 
@@ -360,6 +367,11 @@ def orthogonal_array(k,n,t=2,check=True,availability=False):
       to build such an array. This should be much faster than actually building
       it.
 
+    - ``who_asked`` (internal use only) -- because of the equivalence between
+      OA/TD/MOLS, each of the three constructors calls the others. We must keep
+      track of who calls who in order to avoid infinite loops. ``who_asked`` is
+      the tuple of the other functions that were called before this one.
+
     For more information on orthogonal arrays, see
     :wikipedia:`Orthogonal_array`.
 
@@ -446,12 +458,20 @@ def orthogonal_array(k,n,t=2,check=True,availability=False):
                 l.append(i%n)
             OA = M
 
+    elif (t == 2 and transversal_design not in who_asked and
+          transversal_design(k,n,availability=True, who_asked=who_asked+(orthogonal_array,))):
+        if availability:
+            return True
+        TD = transversal_design(k,n,check=False,who_asked=who_asked+(orthogonal_array,))
+        OA = [[x%n for x in R] for R in TD]
+
     # Section 6.5.1 from [Stinson2004]
-    elif t == 2 and mutually_orthogonal_latin_squares(n,k-2, availability=True):
+    elif (t == 2 and mutually_orthogonal_latin_squares not in who_asked and
+          mutually_orthogonal_latin_squares(n,k-2, availability=True,who_asked=who_asked+(orthogonal_array,))):
         if availability:
             return True
 
-        mols = mutually_orthogonal_latin_squares(n,k-2)
+        mols = mutually_orthogonal_latin_squares(n,k-2,who_asked=who_asked+(orthogonal_array,))
         OA = [[i,j]+[m[i,j] for m in mols]
               for i in range(n) for j in range(n)]