trac #16802: Review

src/sage/combinat/designs/database.py

@@ -3498,10 +3498,16 @@ def OA_33_993():
              [(0,0),(1,0),(3,9),(4,8),(6,1),(9,5)]]},
 (141, 5, 1):
   {(141,): [[0,33,60,92,97],[0,3,45,88,110],[0,18,39,68,139],[0,12,67,75,113],
-               [0,1,15,84,94],[0,7,11,24,30],[0,36,90,116,125]]},
+            [0,1,15,84,94],[0,7,11,24,30],[0,36,90,116,125]]},
 (161, 5, 1):
   {(161,): [[0,19,34,73,80],[0,16,44,71,79],[0,12,33,74,78],[0,13,30,72,77],
             [0,11,36,67,76],[0,18,32,69,75],[0,10,48,68,70],[0,3,29,52,53]]},
+(175, 7, 1):
+  {(7,5,5): [[(0,0,0),(1,0,0),(2,0,0),(3,0,0),(4,0,0),(5,0,0),(6,0,0)],
+             [(0,0,0),(1,1,3),(1,4,2),(2,2,2),(2,3,3),(4,2,0),(4,3,0)],
+             [(0,0,0),(1,3,4),(1,2,1),(2,2,3),(2,3,2),(4,0,2),(4,0,3)],
+             [(0,0,0),(1,1,2),(1,4,3),(2,1,1),(2,4,4),(4,0,1),(4,0,4)],
+             [(0,0,0),(1,3,1),(1,2,4),(2,4,1),(2,1,4),(4,1,0),(4,4,0)]]},
 (201, 5, 1):
   {(201,): [[0,1,45,98,100],[0,3,32,65,89],[0,4,54,70,75],[0,6,49,69,91],[0,7,58,81,95],
             [0,8,34,72,90],[0,9,36,77,96],[0,10,35,83,94],[0,12,40,79,92],[0,15,46,76,93]]},
@@ -3512,22 +3518,16 @@ def OA_33_993():
   {(221,): [[0,1,24,61,116],[0,3,46,65,113],[0,4,73,89,130],[0,5,77,122,124],
             [0,6,39,50,118],[0,7,66,81,94],[0,8,38,64,139],[0,9,29,80,107],
             [0,10,35,93,135],[0,12,34,52,88],[0,14,31,63,84]]},
-# the following one is lemma 2.2 in Abel "Some new BIBDs with block size 7"
-(259, 7, 1):
+
+(259, 7, 1): # the following one is lemma 2.2 in Abel "Some new BIBDs with block size 7"
   {(7,37): [[(0,0),(1,0),(2,0),(3,0),(4,0),(5,0),(6,0)],
             [(0,0),(0,1),(0,6),(1,4),(2,19),(3,25),(6,26)],
             [(0,0),(0,10),(0,23),(2,3),(4,5),(5,1),(6,28)],
             [(0,0),(0,8),(0,26),(1,13),(3,10),(4,30),(5,21)],
             [(0,0),(0,4),(1,25),(1,34),(2,33),(2,35),(4,10)],
             [(0,0),(0,3),(1,26),(2,7),(2,28),(4,17),(4,34)],
             [(0,0),(0,30),(1,7),(1,22),(2,1),(4,21),(4,33)]]},
-(175, 7, 1):
-  {(7,5,5): [[(0,0,0),(1,0,0),(2,0,0),(3,0,0),(4,0,0),(5,0,0),(6,0,0)],
-             [(0,0,0),(1,1,3),(1,4,2),(2,2,2),(2,3,3),(4,2,0),(4,3,0)],
-             [(0,0,0),(1,3,4),(1,2,1),(2,2,3),(2,3,2),(4,0,2),(4,0,3)],
-             [(0,0,0),(1,1,2),(1,4,3),(2,1,1),(2,4,4),(4,0,1),(4,0,4)],
-             [(0,0,0),(1,3,1),(1,2,4),(2,4,1),(2,1,4),(4,1,0),(4,4,0)]]},
-
+
 ##############
 # lambda = 2 #
 ##############
@@ -3593,11 +3593,11 @@ def OA_33_993():
            [9,19,29,37,43,56,59,81]]},
 (153, 9, 2):
   {(3,3,17): [[(0,0,0),(0,1,0),(0,2,0),(1,0,0),(1,1,0),(1,2,0),(2,0,0),(2,1,0),(2,2,0)],
-             [(0,0,0),(0,1,0),(0,2,0),(1,0,0),(1,1,0),(1,2,0),(2,0,0),(2,1,0),(2,2,0)],
-             [(0,0,0),(0,1,1),(0,1,16),(0,2,4),(0,2,13),(1,0,3),(1,0,14),(2,0,5),(2,0,12)],
-             [(0,0,0),(0,1,2),(0,1,15),(0,2,8),(0,2,9),(1,0,6),(1,0,11),(2,0,10),(2,0,7)],
-             [(0,0,0),(0,1,3),(0,1,14),(0,2,12),(0,2,5),(1,0,9),(1,0,8),(2,0,15),(2,0,2)],
-             [(0,0,0),(0,1,6),(0,1,11),(0,2,7),(0,2,10),(1,0,1),(1,0,16),(2,0,13),(2,0,4)]]},
+              [(0,0,0),(0,1,0),(0,2,0),(1,0,0),(1,1,0),(1,2,0),(2,0,0),(2,1,0),(2,2,0)],
+              [(0,0,0),(0,1,1),(0,1,16),(0,2,4),(0,2,13),(1,0,3),(1,0,14),(2,0,5),(2,0,12)],
+              [(0,0,0),(0,1,2),(0,1,15),(0,2,8),(0,2,9),(1,0,6),(1,0,11),(2,0,10),(2,0,7)],
+              [(0,0,0),(0,1,3),(0,1,14),(0,2,12),(0,2,5),(1,0,9),(1,0,8),(2,0,15),(2,0,2)],
+              [(0,0,0),(0,1,6),(0,1,11),(0,2,7),(0,2,10),(1,0,1),(1,0,16),(2,0,13),(2,0,4)]]},
 (181,10, 2):
   {(181,): [[1,7,40,42,51,59,113,125,135,151],
             [19,22,31,35,36,64,74,133,154,156],

src/sage/combinat/designs/difference_family.py

@@ -67,7 +67,16 @@ def group_law(G):
 
 def block_stabilizer(G, B):
     r"""
-    Compute the stabilizer of the block ``B`` in the group ``G``.
+    Compute the left stabilizer of the block ``B`` under the action of ``G``.
+
+    This function return the list of all `x\in G` such that `x\cdot B=B` (as a
+    set).
+
+    INPUT:
+
+    - ``G`` -- a group (additive or multiplicative).
+
+    - ``B`` -- a subset of ``G``.
 
     EXAMPLES::
 
@@ -87,23 +96,25 @@ def block_stabilizer(G, B):
         sage: block_stabilizer(Zmod(45),b)
         [0]
     """
+    if not B:
+        return list(G)
     identity, op, inv = group_law(G)
     b0 = inv(B[0])
     S = []
     for b in B:
         # fun: if we replace +(-b) with -b it completely fails!!
-        bb0 = op(b,b0)
+        bb0 = op(b,b0) # bb0 = b-B[0]
         if all(op(bb0,c) in B for c in B):
             S.append(bb0)
     return S
 
 def is_difference_family(G, D, v=None, k=None, l=None, verbose=False):
     r"""
-    Check wether ``D`` forms a difference family in the additive Abelian group ``G``.
+    Check wether ``D`` forms a difference family in the group ``G``.
 
     INPUT:
 
-    - ``G`` - Additive Abelian group of cardinality ``v``
+    - ``G`` - group of cardinality ``v``
 
     - ``D`` - a set of ``k``-subsets of ``G``
 
@@ -178,77 +189,84 @@ def is_difference_family(G, D, v=None, k=None, l=None, verbose=False):
 
     Glist = list(G)
 
+    D = [map(G,d) for d in D]
+
+    # Check v (and define it if needed)
     if v is None:
         v = len(Glist)
     else:
-        v = int(v)
         if len(Glist) != v:
             if verbose:
-                print "G must have cardinality v (=%d)"%v
+                print "G must have cardinality v (=%d)"%int(v)
             return False
 
+    # Check k (and define it if needed)
     if k is None:
         k = len(D[0])
     else:
         k = int(k)
 
-    b = len(D)
-
     for d in D:
         if len(d) != k:
             if verbose:
                 print "the block {} does not have length {}".format(d,k)
             return False
 
+    # Check l (and define it if needed)
+    #
+    # - nb_diff: the number of pairs (with multiplicity) covered by the BIBD
+    #            generated by the DF.
+    #
+    # - stab: the stabilizer of each set.
     nb_diff = 0
-    stab_sizes = []
+    stab = []
     for d in D:
-        s = len(block_stabilizer(G,map(G,d)))
-        stab_sizes.append(s)
-        nb_diff += k*(k-1) / s
+        s = block_stabilizer(G,d)
+        stab.append(s)
+        nb_diff += k*(k-1) / len(s)
     if l is None:
         if nb_diff % (v-1) != 0:
             if verbose:
-                print "sum(1/s_i) k*(k-1) = {} is not a multiple of (v-1) = {} (stabilizer sizes {})".format(
-                        nb_diff, v-1, stab_sizes)
+                print "sum_i (1/s_i) k*(k-1) = {} is not a multiple of (v-1) = {} (stabilizer sizes {})".format(
+                        nb_diff, v-1, map(len,stab))
             return False
         l = nb_diff // (v-1)
     else:
         if nb_diff != l*(v-1):
             if verbose:
-                print ("the relation sum(1/s_i) *k*(k-1) == l*(v-1) is not "
+                print ("the relation sum_i (1/s_i) *k*(k-1) == l*(v-1) is not "
                 "satisfied (where the s_i are the cardinality of the stabilizer "
                 "for each blocks)")
             return False
 
-    # now we check that every non-identity element of G occurs exactly l-time
-    # as a difference
+    # Check that every x \in G-{0},occurs exactly l times as a difference
     counter = {g: 0 for g in Glist}
     where   = {g: set() for g in Glist}
     del counter[identity]
 
     for i,d in enumerate(D):
-        dd = map(G,d)
-        stab = block_stabilizer(G,dd)
-        ccounter = {}
-        for b in dd:
-            for c in dd:
+        tmp_counter = {}
+        for b in d:
+            for c in d:
                 if b == c:
                     continue
-                gg = mul(b,inv(c))
-                if gg not in ccounter:
-                    ccounter[gg] = 0
+                gg = mul(b,inv(c)) # = b-c
+                if gg not in tmp_counter:
+                    tmp_counter[gg] = 0
                 where[gg].add(i)
-                ccounter[gg] += 1
+                tmp_counter[gg] += 1
 
-        if sum(ccounter.itervalues()) != k*(k-1):
+        if sum(tmp_counter.itervalues()) != k*(k-1):
             if verbose:
                 print "repeated element in the {}-th block {}".format(i,dd)
             return False
 
-        for gg in ccounter:
-            counter[gg] += ccounter[gg]//len(stab)
+        # Normalized number of occurrences added to counter
+        stabi = len(stab[i])
+        for gg in tmp_counter:
+            counter[gg] += tmp_counter[gg]//stabi
 
+    # Check the counter and report any error
     too_few  = []
     too_much = []
     for g in Glist:
@@ -546,23 +564,20 @@ def difference_family(v, k, l=1, existence=False, check=True):
 
         vv, blocks = DF[v,k,l].iteritems().next()
 
+        # Build the group
         from sage.rings.finite_rings.integer_mod_ring import Zmod
         if len(vv) == 1:
-            A = Zmod(vv[0])
+            G = Zmod(vv[0])
         else:
             from sage.categories.cartesian_product import cartesian_product
-            A = cartesian_product([Zmod(i) for i in vv])
+            G = cartesian_product([Zmod(i) for i in vv])
 
-        df = []
-        for b in blocks:
-            b = [A(i) for i in b]
-            df.append(b)
+        df = [[G(i) for i in b] for b in blocks]
 
-        if check and not is_difference_family(A, df, v=v, k=k, l=l):
-            is_difference_family(A,df,v=v,k=k,l=l,verbose=True)
-            raise ValueError("ERROR: pb with ({},{},{})".format(v,k,l))
+        if check:
+            assert is_difference_family(G, df, v=v, k=k, l=l), "Sage built an invalid ({},{},{})-DF!".format(v,k,l)
 
-        return A,df
+        return G,df
 
     e = k*(k-1)
     t = l*(v-1) // e  # number of blocks