Removed the use of baseswap from subgroup search; minor changes.

Due to a bug in the version of subgroup_search using baseswap
to perform the base changes, the way basic stabilizers are
obtained is now via the stabilizer() function. The function
_insert_point_in_base was remove from sympy.combinatorics.util
since it's not needed right now.

Also, made the incremental Schreier-Sims algorithm exclude the
identity element as a generator and wrote some more tests
for subgroup_search

sympy/combinatorics/perm_groups.py

@@ -8,8 +8,7 @@
 from sympy.ntheory import isprime, sieve
 from sympy.combinatorics.util import _check_cycles_alt_sym,\
 _distribute_gens_by_base, _orbits_transversals_from_bsgs,\
-_handle_precomputed_bsgs, _base_ordering, _strong_gens_from_distr, _strip,\
-_insert_point_in_base
+_handle_precomputed_bsgs, _base_ordering, _strong_gens_from_distr, _strip
 
 def _smallest_change(h, alpha):
     """
@@ -2086,12 +2085,15 @@ def schreier_sims_incremental(self, base=None, gens=None):
         if base is None:
             base = []
         if gens is None:
-            gens = self.generators
+            gens = self.generators[:]
         base_len = len(base)
         degree = self.degree
+        identity = _new_from_array_form(range(degree))
         # handle the trivial group
-        if gens == [_new_from_array_form(range(degree))]:
+        if gens == [identity]:
             return base, gens
+        # remove the identity as a generator
+        gens = [x for x in gens if x != identity]
         # make sure no generator fixes all base points
         for gen in gens:
             if [gen(x) for x in base] == [x for x in base]:
@@ -2506,6 +2508,8 @@ def subgroup_search(self, prop, base=None, strong_gens=None, tests=None, init_su
         # initialize computed words
         computed_words = [identity]*base_len
         # line 8: main loop
+        # stabilizers
+        stabilizers_res = [PermutationGroup(gens) for gens in res_distr_gens]
         while True:
             # apply all the tests
             while l < base_len - 1 and\
@@ -2515,10 +2519,9 @@ def subgroup_search(self, prop, base=None, strong_gens=None, tests=None, init_su
                   tests[l](computed_words[base_len - 1]):
                 # line 11: change the (partial) base of K
                 new_point = computed_words[l](base[l])
-                # this function maintains a partial BSGS structure up to position l
-                _insert_point_in_base(res, res_base, res_strong_gens, l, new_point, distr_gens=res_distr_gens, basic_orbits=res_basic_orbits, transversals=res_transversals)
-                # find the l+1-th basic stabilizer
-                new_stab = PermutationGroup(res_distr_gens[l + 1])
+                res_base[l] = new_point
+                stabilizers_res[l + 1] = stabilizers_res[l].stabilizer(new_point)
+                new_stab = stabilizers_res[l + 1]
                 # line 12: calculate minimal orbit representatives for the l+1-th basic stabilizer
                 orbits = new_stab.orbits()
                 reps = []
@@ -2573,7 +2576,8 @@ def subgroup_search(self, prop, base=None, strong_gens=None, tests=None, init_su
                 res_basic_orbits_init_base = res_basic_orbits[:]
                 res_strong_gens = res_strong_gens_init[:]
                 # line 21: recalculate orbit representatives
-                stab_f = PermutationGroup(res_distr_gens[f])
+                stabilizers_res = [PermutationGroup(gens) for gens in res_distr_gens]
+                stab_f = stabilizers_res[f]
                 temp_orbits = stab_f.orbits()
                 reps = []
                 for orbit in orbits:
@@ -2594,7 +2598,7 @@ def subgroup_search(self, prop, base=None, strong_gens=None, tests=None, init_su
                 f = l
                 c[l] = 0
                 # line 27
-                stab_f = PermutationGroup(res_distr_gens[f])
+                stab_f = stabilizers_res[f]
                 temp_orbits = stab_f.orbits()
                 reps = []
                 for orbit in orbits:

sympy/combinatorics/tests/test_perm_groups.py

@@ -405,11 +405,29 @@ def test_schreier_sims_incremental():
 def test_subgroup_search():
     prop_true = lambda x: True
     prop_fix_points = lambda x: [x(point) for point in points] == points
-    prop_commutes_with_g = lambda x: x*g == g*x
+    prop_comm_g = lambda x: x*g == g*x
     prop_even = lambda x: x.is_even
-    S = SymmetricGroup(11)
-    A = AlternatingGroup(11)
-    C = CyclicGroup(11)
-    assert S.subgroup_search(prop_true) == S
-    assert S.subgroup_search(prop_true, init_subgroup = C) == S
-    assert S.subgroup_search(prop_even) == A
+    for i in range(10, 17, 2):
+        S = SymmetricGroup(i)
+        A = AlternatingGroup(i)
+        C = CyclicGroup(i)
+        Sym = S.subgroup_search(prop_true)
+        assert Sym == S
+        Alt = S.subgroup_search(prop_even)
+        assert Alt == A
+        Sym = S.subgroup_search(prop_true, init_subgroup=C)
+        assert Sym == S
+        points = [7]
+        assert S.stabilizer(7) == S.subgroup_search(prop_fix_points)
+        points = [3, 4]
+        assert S.stabilizer(3).stabilizer(4) == S.subgroup_search(prop_fix_points)
+        points = [3, 5]
+        fix35 = A.subgroup_search(prop_fix_points)
+        points = [5]
+        fix5 = A.subgroup_search(prop_fix_points)
+        assert A.subgroup_search(prop_fix_points, init_subgroup=fix35) == fix5
+        base, strong_gens = A.schreier_sims_incremental()
+        g = A.generators[0]
+        comm_g = A.subgroup_search(prop_comm_g, base=base, strong_gens=strong_gens)
+        assert _verify_bsgs(comm_g, base, comm_g.generators) == True
+        assert [prop_comm_g(gen) == True for gen in comm_g.generators]

sympy/combinatorics/util.py

@@ -236,57 +236,6 @@ def _handle_precomputed_bsgs(base, strong_gens, transversals=None,\
                 basic_orbits[i] = transversals[i].keys()
     return transversals, basic_orbits, distr_gens
 
-def _insert_point_in_base(group, base, strong_gens, pos, point, distr_gens=None, basic_orbits=None, transversals=None):
-    from sympy.combinatorics.perm_groups import PermutationGroup
-    # initialize basic group properties and BSGS structures
-    base_len = len(base)
-    degree = group.degree
-    identity = _new_from_array_form(range(degree))
-    transversals, basic_orbits, distr_gens = _handle_precomputed_bsgs(base, strong_gens, transversals=transversals, basic_orbits=basic_orbits, distr_gens=distr_gens)
-    # cut the base at position pos and append the new point
-    partial_base = base[: pos + 1]
-    partial_base.append(point)
-    # cut the generators for the stabilizer chain and amend them accordingly
-    if pos == base_len - 1:
-        partial_distr_gens = distr_gens[: pos + 1]
-        partial_distr_gens.append([identity])
-    else:
-        partial_distr_gens = distr_gens[: pos + 2]
-    # cut the basic orbits and transversals and amend them accordingly
-    partial_basic_orbits = basic_orbits[: pos + 1]
-    partial_transversals = transversals[: pos + 1]
-    last_stab = PermutationGroup(partial_distr_gens[pos + 1])
-    last_transversal = dict(last_stab.orbit_transversal(point, pairs=True))
-    last_orbit = last_transversal.keys()
-    partial_basic_orbits.append(last_orbit)
-    partial_transversals.append(last_transversal)
-    # baseswap with the partial BSGS structures. Notice that we need only
-    # the orbit and transversal of the new point under the last stabilizer
-    new_base, new_strong_gens = group.baseswap(partial_base, strong_gens, pos, randomized=False, transversals=partial_transversals, basic_orbits=partial_basic_orbits, distr_gens=partial_distr_gens)
-    # amend the basic orbits and transversals
-    stab_pos = PermutationGroup(distr_gens[pos])
-    new_transversal = dict(stab_pos.orbit_transversal(point, pairs=True))
-    transversals[pos] = new_transversal
-    basic_orbits[pos] = new_transversal.keys()
-    # amend the distributed generators if necessary
-    if pos != base_len - 1:
-        new_stab_gens = []
-        for gen in new_strong_gens:
-            if [gen(point) for point in new_base[: pos + 1]] == [point for point in new_base[: pos + 1]]:
-                new_stab_gens.append(gen)
-        distr_gens[pos + 1] = new_stab_gens
-    # return the new partial base and partial strong generating set
-    new_base.pop()
-    new_base = new_base + base[pos + 1 :]
-    while len(base) != 0:
-        base.pop()
-    for point in new_base:
-        base.append(point)
-    while len(strong_gens) != 0:
-        strong_gens.pop()
-    for gen in new_strong_gens:
-        strong_gens.append(gen)
-
 def _orbits_transversals_from_bsgs(base, distr_gens,\
                                    transversals_only=False):
     """