Skip to content
This repository

Incremental Schreier-Sims algorithm, and subgroup searching. #1454

Merged
merged 13 commits into from over 1 year ago

5 participants

Aleksandar Makelov Don't Add Me To Your Organization a.k.a The Travis Bot Stefan Krastanov david joyner Aaron Meurer
Aleksandar Makelov

Implemented a version of the Schreier-Sims algorithm that takes a sequence of points and a generating set for a group and extends them to a base and strong generating set.

Implemented a procedure to remove redundant generators from a strong generating set.

Implemented the procedure subgroup_search which finds generators for the subgroup of all elements satisfying a given property in a larger group.

Don't Add Me To Your Organization a.k.a The Travis Bot

This pull request passes (merged a6253a17 into 02daa7d).

Stefan Krastanov
Collaborator

SymPy Bot Summary: :red_circle: There were test failures.

@amakelov: Please fix the test failures.

Test command: setup.py test
master hash: 02daa7d
branch hash: a6253a1

Interpreter 1: :red_circle: There were test failures.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYnbYiDA

Interpreter 2: :red_circle: There were test failures.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY3PQiDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYmI8iDA

Build HTML Docs: :red_circle: There were test failures.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYn90iDA

Automatic review by SymPy Bot.

Don't Add Me To Your Organization a.k.a The Travis Bot

This pull request passes (merged 02198cca into 02daa7d).

Stefan Krastanov
Collaborator

SymPy Bot Summary: :red_circle: There were test failures.

@amakelov: Please fix the test failures.

Test command: setup.py test
master hash: bff3595
branch hash: 02198cc

Interpreter 1: :red_circle: There were test failures.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY8MUiDA

Interpreter 2: :red_circle: There were test failures.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYiIwjDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY4fQiDA

Build HTML Docs: :red_circle: There were test failures.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY1M0iDA

Automatic review by SymPy Bot.

added some commits July 29, 2012
Aleksandar Makelov Incremental Schreier-Sims algorithm.
It is used to extend a sequence of points and a generating set
to a base and strong generating set relative to it.
c592584
Aleksandar Makelov Tests for incremental Schreier-Sims algorithm. 2090a19
Aleksandar Makelov Remove redundant strong generators from a strong generating set.
Via the function _remove_gens in sympy/combinatorics/util.py
b852d89
Aleksandar Makelov Subgroup search and tests for _remove_gens.
The function subgroup_search is used to find a strong generating
set for the subgroup of all elements satisfying a given property
within a group.
454268c
Aleksandar Makelov Added tests for subgroup search. be08509
Aleksandar Makelov 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
0cdee22
Aleksandar Makelov Removed unnecessary code from subgroup_search. 1ad8444
Aleksandar Makelov Added docstrings for the new functions. fabddb8
Aleksandar Makelov Minor style improvements. 8820277
Aleksandar Makelov Got rid of lines longer than 80 characters. ea648ae
Aleksandar Makelov Renamed all variables `distr_gens` to `strong_gens_distr`.
For more clarity, as suggested by my GSoC mentor David. This
variable is usually used to denote a list of strong generators
distributed by membership in basic stabilizers.
03f2c95
Aleksandar Makelov Fixed some issues with the docs. b4cc733
Aleksandar Makelov Removed side effects from schreier_sims_incremental. 4dffbde
Don't Add Me To Your Organization a.k.a The Travis Bot

This pull request passes (merged 4dffbde into 65b6582).

Stefan Krastanov
Collaborator

SymPy Bot Summary: :eight_spoked_asterisk: All tests have passed.

Test command: setup.py test
master hash: 3c2c3c7
branch hash: 4dffbde

Interpreter 1: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYhsYiDA

Interpreter 2: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY9JsjDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY8f8hDA

Build HTML Docs: :eight_spoked_asterisk: All tests have passed.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYp4wjDA

Automatic review by SymPy Bot.

david joyner

This patch looks good to me and it passes all tests.

Aaron Meurer
Owner

OK, I'll merge it then.

Aaron Meurer asmeurer merged commit c7a4a4a into from August 07, 2012
Aaron Meurer asmeurer closed this August 07, 2012
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Showing 13 unique commits by 1 author.

Aug 06, 2012
Aleksandar Makelov Incremental Schreier-Sims algorithm.
It is used to extend a sequence of points and a generating set
to a base and strong generating set relative to it.
c592584
Aleksandar Makelov Tests for incremental Schreier-Sims algorithm. 2090a19
Aleksandar Makelov Remove redundant strong generators from a strong generating set.
Via the function _remove_gens in sympy/combinatorics/util.py
b852d89
Aleksandar Makelov Subgroup search and tests for _remove_gens.
The function subgroup_search is used to find a strong generating
set for the subgroup of all elements satisfying a given property
within a group.
454268c
Aleksandar Makelov Added tests for subgroup search. be08509
Aleksandar Makelov 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
0cdee22
Aleksandar Makelov Removed unnecessary code from subgroup_search. 1ad8444
Aleksandar Makelov Added docstrings for the new functions. fabddb8
Aleksandar Makelov Minor style improvements. 8820277
Aleksandar Makelov Got rid of lines longer than 80 characters. ea648ae
Aleksandar Makelov Renamed all variables `distr_gens` to `strong_gens_distr`.
For more clarity, as suggested by my GSoC mentor David. This
variable is usually used to denote a list of strong generators
distributed by membership in basic stabilizers.
03f2c95
Aleksandar Makelov Fixed some issues with the docs. b4cc733
Aleksandar Makelov Removed side effects from schreier_sims_incremental. 4dffbde
This page is out of date. Refresh to see the latest.
1  doc/src/modules/combinatorics/index.rst
Source Rendered
@@ -16,3 +16,4 @@ Contents
16 16
    subsets.rst
17 17
    graycode.rst
18 18
    named_groups.rst
  19
+   util.rst
12  doc/src/modules/combinatorics/util.txt → doc/src/modules/combinatorics/util.rst
Source Rendered
@@ -5,18 +5,20 @@ Utilities
5 5
 
6 6
 .. module:: sympy.combinatorics.util
7 7
 
8  
-.. autofunction:: _check_cycles_alt_sym
  8
+.. autofunction:: _base_ordering
9 9
 
10  
-.. autofunction:: _strip
  10
+.. autofunction:: _check_cycles_alt_sym
11 11
 
12 12
 .. autofunction:: _distribute_gens_by_base
13 13
 
14  
-.. autofunction:: _strong_gens_from_distr
  14
+.. autofunction:: _handle_precomputed_bsgs
15 15
 
16 16
 .. autofunction:: _orbits_transversals_from_bsgs
17 17
 
18  
-.. autofunction:: _handle_precomputed_bsgs
  18
+.. autofunction:: _remove_gens
19 19
 
20  
-.. autofunction:: _base_ordering
  20
+.. autofunction:: _strip
  21
+
  22
+.. autofunction:: _strong_gens_from_distr
21 23
 
22 24
 .. autofunction:: _verify_bsgs
4  sympy/combinatorics/named_groups.py
@@ -70,13 +70,13 @@ def AlternatingGroup(n):
70 70
 
71 71
     """
72 72
     # small cases are special
73  
-    if n == 1 or n == 2:
  73
+    if n in (1, 2):
74 74
         return PermutationGroup([Permutation([0])])
75 75
 
76 76
     a = range(n)
77 77
     a[0], a[1], a[2] = a[1], a[2], a[0]
78 78
     gen1 = _new_from_array_form(a)
79  
-    if n % 2 == 1:
  79
+    if n % 2:
80 80
         a = range(1, n)
81 81
         a.append(0)
82 82
         gen2 = _new_from_array_form(a)
441  sympy/combinatorics/perm_groups.py
@@ -475,7 +475,7 @@ def _random_pr_init(self, r, n, _random_prec_n=None):
475 475
         Notes
476 476
         =====
477 477
 
478  
-        THIS FUNCTION HAS SIDE EFFECTS: it changes the attribute
  478
+        XXX THIS FUNCTION HAS SIDE EFFECTS: it changes the attribute
479 479
         self._random_gens
480 480
 
481 481
         See Also
@@ -619,8 +619,8 @@ def base(self):
619 619
             self.schreier_sims()
620 620
         return self._base
621 621
 
622  
-    def baseswap(self, base, strong_gens, pos, randomized=True,\
623  
-                 transversals=None, basic_orbits=None, distr_gens=None):
  622
+    def baseswap(self, base, strong_gens, pos, randomized=False,\
  623
+                 transversals=None, basic_orbits=None, strong_gens_distr=None):
624 624
         r"""
625 625
         Swap two consecutive base points in a base and strong generating set.
626 626
 
@@ -639,7 +639,7 @@ def baseswap(self, base, strong_gens, pos, randomized=True,\
639 639
         ``randomized`` - switch between randomized and deterministic version
640 640
         ``transversals`` - transversals for the basic orbits, if known
641 641
         ``basic_orbits`` - basic orbits, if known
642  
-        ``distr_gens`` - strong generators distributed by basic stabilizers,
  642
+        ``strong_gens_distr`` - strong generators distributed by basic stabilizers,
643 643
         if known
644 644
 
645 645
         Returns
@@ -682,12 +682,12 @@ def baseswap(self, base, strong_gens, pos, randomized=True,\
682 682
         """
683 683
         # construct the basic orbits, generators for the stabilizer chain
684 684
         # and transversal elements from whatever was provided
685  
-        transversals, basic_orbits, distr_gens =\
  685
+        transversals, basic_orbits, strong_gens_distr =\
686 686
         _handle_precomputed_bsgs(base, strong_gens, transversals,\
687  
-                                 basic_orbits, distr_gens)
  687
+                                 basic_orbits, strong_gens_distr)
688 688
         base_len = len(base)
689 689
         degree = self.degree
690  
-        stab_pos = PermutationGroup(distr_gens[pos])
  690
+        stab_pos = PermutationGroup(strong_gens_distr[pos])
691 691
         # size of orbit of base[pos] under the stabilizer we seek to insert
692 692
         # in the stabilizer chain at position pos + 1
693 693
         size = len(basic_orbits[pos])*len(basic_orbits[pos + 1])\
@@ -696,7 +696,7 @@ def baseswap(self, base, strong_gens, pos, randomized=True,\
696 696
         if pos + 2 > base_len - 1:
697 697
             T = []
698 698
         else:
699  
-            T = distr_gens[pos + 2][:]
  699
+            T = strong_gens_distr[pos + 2][:]
700 700
         if T == []:
701 701
             current_group = PermGroup([_new_from_array_form(range(degree))])
702 702
         else:
@@ -733,7 +733,7 @@ def baseswap(self, base, strong_gens, pos, randomized=True,\
733 733
                         current_group = PermutationGroup(T)
734 734
                         Gamma = Gamma - current_group.orbit(base[pos])
735 735
         # build the new base and strong generating set
736  
-        strong_gens_new_distr = distr_gens[:]
  736
+        strong_gens_new_distr = strong_gens_distr[:]
737 737
         strong_gens_new_distr[pos + 1] = T
738 738
         base_new = base[:]
739 739
         base_new[pos], base_new[pos + 1] = base_new[pos + 1], base_new[pos]
@@ -803,9 +803,9 @@ def basic_stabilizers(self):
803 803
             self.schreier_sims()
804 804
         strong_gens = self._strong_gens
805 805
         base = self._base
806  
-        distr_gens = _distribute_gens_by_base(base, strong_gens)
  806
+        strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
807 807
         basic_stabilizers = []
808  
-        for gens in distr_gens:
  808
+        for gens in strong_gens_distr:
809 809
             basic_stabilizers.append(PermutationGroup(gens))
810 810
         return basic_stabilizers
811 811
 
@@ -1732,7 +1732,7 @@ def orbit(self, alpha, action='tuples'):
1732 1732
                         orb.append(temp)
1733 1733
                         used[temp] = True
1734 1734
             return set(orb)
1735  
-        if action == 'tuples':
  1735
+        elif action == 'tuples':
1736 1736
             alpha = tuple(alpha)
1737 1737
             orb = [alpha]
1738 1738
             used = set([alpha])
@@ -1744,7 +1744,7 @@ def orbit(self, alpha, action='tuples'):
1744 1744
                         orb.append(temp)
1745 1745
                         used.add(temp)
1746 1746
             return set(orb)
1747  
-        if action == 'sets':
  1747
+        elif action == 'sets':
1748 1748
             alpha = frozenset(alpha)
1749 1749
             orb = [alpha]
1750 1750
             used = set([alpha])
@@ -2081,6 +2081,150 @@ def schreier_sims(self):
2081 2081
         self._transversals = transversals
2082 2082
         self._basic_orbits = basic_orbits
2083 2083
 
  2084
+    def schreier_sims_incremental(self, base=None, gens=None):
  2085
+        """
  2086
+        Extend a sequence of points and generating set to a base and strong
  2087
+        generating set.
  2088
+
  2089
+        Parameters
  2090
+        ==========
  2091
+
  2092
+        ``base`` - the sequence of points to be extended to a base. Optional
  2093
+        parameter with default value ``[]``
  2094
+        ``gens`` - generating set to be extended to a strong generating set
  2095
+        relative to the base obtained. Optional parameter with default value
  2096
+        ``self.generators``
  2097
+
  2098
+        Returns
  2099
+        =======
  2100
+
  2101
+        ``(base, strong_gens)`` where ``base`` is the base obtained, and
  2102
+        ``strong_gens`` is the strong generating set relative to it. The
  2103
+        original parameters ``base``, ``gens`` remain unchanged.
  2104
+
  2105
+        Examples
  2106
+        ========
  2107
+
  2108
+        >>> from sympy.combinatorics.named_groups import AlternatingGroup
  2109
+        >>> from sympy.combinatorics.perm_groups import PermutationGroup
  2110
+        >>> from sympy.combinatorics.util import _verify_bsgs
  2111
+        >>> A = AlternatingGroup(7)
  2112
+        >>> base = [2, 3]
  2113
+        >>> seq = [2, 3]
  2114
+        >>> base, strong_gens = A.schreier_sims_incremental(base=seq)
  2115
+        >>> _verify_bsgs(A, base, strong_gens)
  2116
+        True
  2117
+        >>> base[:2]
  2118
+        [2, 3]
  2119
+
  2120
+        Notes
  2121
+        =====
  2122
+
  2123
+        This version of the Schreier-Sims algorithm runs in polynomial time.
  2124
+        There are certain assumptions in the implementation - if the trivial
  2125
+        group is provided, ``base`` and ``gens`` are returned immediately,
  2126
+        as any sequence of points is a base for the trivial group. If the
  2127
+        identity is present in the generators ``gens``, it is removed as
  2128
+        it is a redundant generator.
  2129
+        The implementation is described in [1],pp.90-93.
  2130
+
  2131
+        See Also
  2132
+        ========
  2133
+
  2134
+        schreier_sims, schreier_sims_random
  2135
+
  2136
+        """
  2137
+        if base is None:
  2138
+            base = []
  2139
+        if gens is None:
  2140
+            gens = self.generators[:]
  2141
+        base_len = len(base)
  2142
+        degree = self.degree
  2143
+        identity = _new_from_array_form(range(degree))
  2144
+        # handle the trivial group
  2145
+        if gens == [identity]:
  2146
+            return base, gens
  2147
+        # prevent side effects
  2148
+        _base, _gens = base[:], gens[:]
  2149
+        # remove the identity as a generator
  2150
+        _gens = [x for x in _gens if x != identity]
  2151
+        # make sure no generator fixes all base points
  2152
+        for gen in _gens:
  2153
+            if [gen(x) for x in _base] == [x for x in _base]:
  2154
+                new = 0
  2155
+                while gen(new) == new:
  2156
+                    new += 1
  2157
+                _base.append(new)
  2158
+                base_len += 1
  2159
+        # distribute generators according to basic stabilizers
  2160
+        strong_gens_distr = _distribute_gens_by_base(_base, _gens)
  2161
+        # initialize the basic stabilizers, basic orbits and basic transversals
  2162
+        stabs = {}
  2163
+        orbs = {}
  2164
+        transversals = {}
  2165
+        for i in xrange(base_len):
  2166
+            stabs[i] = PermutationGroup(strong_gens_distr[i])
  2167
+            transversals[i] = dict(stabs[i].orbit_transversal(_base[i],\
  2168
+                                                              pairs=True))
  2169
+            orbs[i] = transversals[i].keys()
  2170
+        # main loop: amend the stabilizer chain until we have generators
  2171
+        # for all stabilizers
  2172
+        i = base_len - 1
  2173
+        while i >= 0:
  2174
+            # this flag is used to continue with the main loop from inside
  2175
+            # a nested loop
  2176
+            continue_i = False
  2177
+            # test the generators for being a strong generating set
  2178
+            for beta in orbs[i]:
  2179
+                u_beta = transversals[i][beta]
  2180
+                for gen in strong_gens_distr[i]:
  2181
+                    u_beta_gen = transversals[i][gen(beta)]
  2182
+                    if gen*u_beta != u_beta_gen:
  2183
+                        # test if the schreier generator is in the i+1-th
  2184
+                        # would-be basic stabilizer
  2185
+                        y = True
  2186
+                        schreier_gen = (~u_beta_gen)*gen*u_beta
  2187
+                        h, j = _strip(schreier_gen, _base, orbs, transversals)
  2188
+                        if j <= base_len:
  2189
+                            # new strong generator h at level j
  2190
+                            y = False
  2191
+                        elif h != _new_from_array_form(range(degree)):
  2192
+                            # h fixes all base points
  2193
+                            y = False
  2194
+                            moved = 0
  2195
+                            while h(moved) == moved:
  2196
+                                moved += 1
  2197
+                            _base.append(moved)
  2198
+                            base_len += 1
  2199
+                            strong_gens_distr.append([])
  2200
+                        if y == False:
  2201
+                            # if a new strong generator is found, update the
  2202
+                            # data structures and start over
  2203
+                            for l in range(i + 1, j):
  2204
+                                strong_gens_distr[l].append(h)
  2205
+                                stabs[l] = PermutationGroup(strong_gens_distr[l])
  2206
+                                transversals[l] =\
  2207
+                                dict(stabs[l].orbit_transversal(_base[l],\
  2208
+                                                                pairs=True))
  2209
+                                orbs[l] = transversals[l].keys()
  2210
+                            i = j - 1
  2211
+                            # continue main loop using the flag
  2212
+                            continue_i = True
  2213
+                    if continue_i == True:
  2214
+                        break
  2215
+                if continue_i == True:
  2216
+                    break
  2217
+            if continue_i == True:
  2218
+                continue
  2219
+            i -= 1
  2220
+        # build the strong generating set
  2221
+        strong_gens = []
  2222
+        for gens in strong_gens_distr:
  2223
+            for gen in gens:
  2224
+                if gen not in strong_gens:
  2225
+                    strong_gens.append(gen)
  2226
+        return _base, strong_gens
  2227
+
2084 2228
     def schreier_sims_random(self, base=None, gens=None, consec_succ=10,\
2085 2229
                              _random_prec=None):
2086 2230
         r"""
@@ -2158,13 +2302,13 @@ def schreier_sims_random(self, base=None, gens=None, consec_succ=10,\
2158 2302
                 base.append(new)
2159 2303
                 base_len += 1
2160 2304
         # distribute generators according to basic stabilizers
2161  
-        distr_gens = _distribute_gens_by_base(base, gens)
  2305
+        strong_gens_distr = _distribute_gens_by_base(base, gens)
2162 2306
         # initialize the basic stabilizers, basic transversals and basic orbits
2163 2307
         stabs = {}
2164 2308
         transversals = {}
2165 2309
         orbs = {}
2166 2310
         for i in xrange(base_len):
2167  
-            stabs[i] = PermutationGroup(distr_gens[i])
  2311
+            stabs[i] = PermutationGroup(strong_gens_distr[i])
2168 2312
             transversals[i] = dict(stabs[i].orbit_transversal(base[i],\
2169 2313
                                                               pairs=True))
2170 2314
             orbs[i] = transversals[i].keys()
@@ -2189,13 +2333,13 @@ def schreier_sims_random(self, base=None, gens=None, consec_succ=10,\
2189 2333
                     moved += 1
2190 2334
                 base.append(moved)
2191 2335
                 base_len += 1
2192  
-                distr_gens.append([])
  2336
+                strong_gens_distr.append([])
2193 2337
             # if the element doesn't sift, amend the strong generators and
2194 2338
             # associated stabilizers and orbits
2195 2339
             if y == False:
2196 2340
                 for l in range(1, j):
2197  
-                    distr_gens[l].append(h)
2198  
-                    stabs[l] = PermutationGroup(distr_gens[l])
  2341
+                    strong_gens_distr[l].append(h)
  2342
+                    stabs[l] = PermutationGroup(strong_gens_distr[l])
2199 2343
                     transversals[l] = dict(stabs[l].orbit_transversal(base[l],\
2200 2344
                                                                     pairs=True))
2201 2345
                     orbs[l] = transversals[l].keys()
@@ -2203,8 +2347,8 @@ def schreier_sims_random(self, base=None, gens=None, consec_succ=10,\
2203 2347
             else:
2204 2348
                 c += 1
2205 2349
         # build the strong generating set
2206  
-        strong_gens = distr_gens[0][:]
2207  
-        for gen in distr_gens[1]:
  2350
+        strong_gens = strong_gens_distr[0][:]
  2351
+        for gen in strong_gens_distr[1]:
2208 2352
             if gen not in strong_gens:
2209 2353
                 strong_gens.append(gen)
2210 2354
         return base, strong_gens
@@ -2361,6 +2505,263 @@ def strong_gens(self):
2361 2505
             self.schreier_sims()
2362 2506
         return self._strong_gens
2363 2507
 
  2508
+    def subgroup_search(self, prop, base=None, strong_gens=None, tests=None,\
  2509
+                        init_subgroup=None):
  2510
+        """
  2511
+        Find the subgroup of all elements satisfying the property ``prop``.
  2512
+
  2513
+        This is done by a depth-first search with respect to base images that
  2514
+        uses several tests to prune the search tree.
  2515
+
  2516
+        Parameters
  2517
+        ==========
  2518
+
  2519
+        ``prop`` - the property to be used. Has to be callable on group
  2520
+        elements and always return ``True`` or ``False``. It is assumed that
  2521
+        all group elements satisfying ``prop`` indeed form a subgroup.
  2522
+        ``base`` - a base for the supergroup.
  2523
+        ``strong_gens`` - a strong generating set for the supergroup.
  2524
+        ``tests`` - list of callables of length equal to the length of ``base``.
  2525
+        These are used to rule out group elements by partial base images, so
  2526
+        that ``tests[l](g)`` returns False if the element ``g`` is known not
  2527
+        to satisfy prop base on where g sends the first ``l + 1`` base points.
  2528
+        ``init_subgroup`` - if a subgroup of the saught group is known in
  2529
+        advance, it can be passed to the function as this parameter.
  2530
+
  2531
+        Returns
  2532
+        =======
  2533
+
  2534
+        The subgroup of all elements satisfying ``prop``. The generating set
  2535
+        for this group is guaranteed to be a strong generating set relative to
  2536
+        the base ``base``.
  2537
+
  2538
+        Examples
  2539
+        ========
  2540
+
  2541
+        >>> from sympy.combinatorics.named_groups import (SymmetricGroup,
  2542
+        ... AlternatingGroup)
  2543
+        >>> from sympy.combinatorics.perm_groups import PermutationGroup
  2544
+        >>> from sympy.combinatorics.util import _verify_bsgs
  2545
+        >>> S = SymmetricGroup(7)
  2546
+        >>> prop_even = lambda x: x.is_even
  2547
+        >>> base, strong_gens = S.schreier_sims_incremental()
  2548
+        >>> G = S.subgroup_search(prop_even, base=base, strong_gens=strong_gens)
  2549
+        >>> G == AlternatingGroup(7)
  2550
+        True
  2551
+        >>> _verify_bsgs(G, base, G.generators)
  2552
+        True
  2553
+
  2554
+        Notes
  2555
+        =====
  2556
+
  2557
+        This function is extremely lenghty and complicated and will require
  2558
+        some careful attention. The implementation is described in
  2559
+        [1],pp.114-117, and the comments for the code here follow the lines
  2560
+        of the pseudocode in the book for clarity.
  2561
+        The complexity is exponential in general, since the search process by
  2562
+        itself visits all members of the supergroup. However, there are a lot
  2563
+        of tests which are used to prune the search tree, and users can define
  2564
+        their own tests via the ``tests`` parameter, so in practice, and for
  2565
+        some computations, it's not terrible.
  2566
+        A crucial part in the procedure is the frequent base change performed
  2567
+        (this is line 11 in the pseudocode) in order to obtain a new basic
  2568
+        stabilizer. The book mentiones that this can be done by using
  2569
+        ``.baseswap(...)``, however the current imlementation uses a more
  2570
+        straightforward way to find the next basic stabilizer - calling the
  2571
+        function ``.stabilizer(...)`` on the previous basic stabilizer.
  2572
+
  2573
+        """
  2574
+        # initialize BSGS and basic group properties
  2575
+        if base is None:
  2576
+            base, strong_gens = self.schreier_sims_incremental()
  2577
+        base_len = len(base)
  2578
+        degree = self.degree
  2579
+        identity = _new_from_array_form(range(degree))
  2580
+        base_ordering = _base_ordering(base, degree)
  2581
+        # add an element larger than all points
  2582
+        base_ordering.append(degree)
  2583
+        # add an element smaller than all points
  2584
+        base_ordering.append(-1)
  2585
+        # compute BSGS-related structures
  2586
+        strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
  2587
+        basic_orbits, transversals = _orbits_transversals_from_bsgs(base,\
  2588
+                                                                    strong_gens_distr)
  2589
+        # handle subgroup initialization and tests
  2590
+        if init_subgroup is None:
  2591
+            init_subgroup = PermutationGroup([identity])
  2592
+        if tests is None:
  2593
+            trivial_test = lambda x: True
  2594
+            tests = []
  2595
+            for i in xrange(base_len):
  2596
+                tests.append(trivial_test)
  2597
+        # line 1: more initializations.
  2598
+        res = init_subgroup
  2599
+        f = base_len - 1
  2600
+        l = base_len - 1
  2601
+        # line 2: set the base for K to the base for G
  2602
+        res_base = base[:]
  2603
+        # line 3: compute BSGS and related structures for K
  2604
+        res_base, res_strong_gens = res.schreier_sims_incremental(base=res_base)
  2605
+        res_strong_gens_distr = _distribute_gens_by_base(res_base, res_strong_gens)
  2606
+        res_basic_orbits_init_base =\
  2607
+        [PermutationGroup(res_strong_gens_distr[i]).orbit(res_base[i])\
  2608
+         for i in range(base_len)]
  2609
+        # initialize orbit representatives
  2610
+        orbit_reps = [None]*base_len
  2611
+        # line 4: orbit representatives for f-th basic stabilizer of K
  2612
+        stab_f = PermutationGroup(res_strong_gens_distr[f])
  2613
+        orbits = stab_f.orbits()
  2614
+        reps = []
  2615
+        for orbit in orbits:
  2616
+            # get the minimal element in the base ordering
  2617
+            rep = min(orbit, key = lambda point: base_ordering[point])
  2618
+            reps.append(rep)
  2619
+        orbit_reps[f] = reps
  2620
+        # line 5: remove the base point from the representatives to avoid
  2621
+        # getting the identity element as a generator for K
  2622
+        orbit_reps[f].remove(base[f])
  2623
+        # line 6: more initializations
  2624
+        c = [0]*base_len
  2625
+        u = [identity]*base_len
  2626
+        sorted_orbits = [None]*base_len
  2627
+        for i in range(base_len):
  2628
+            sorted_orbits[i] = basic_orbits[i][:]
  2629
+            sorted_orbits[i].sort(key = lambda point: base_ordering[point])
  2630
+        # line 7: initializations
  2631
+        mu = [None]*base_len
  2632
+        nu = [None]*base_len
  2633
+        # this corresponds to the element smaller than all points
  2634
+        mu[l] = degree + 1
  2635
+        temp_index = len(basic_orbits[l])+1-len(res_basic_orbits_init_base[l])
  2636
+        if temp_index >= len(basic_orbits[l]):
  2637
+            # this corresponds to the element larger than all points
  2638
+            nu[l] = base_ordering[degree]
  2639
+        else:
  2640
+            nu[l] = sorted_orbits[l][temp_index]
  2641
+        # initialize computed words
  2642
+        computed_words = [identity]*base_len
  2643
+        # line 8: main loop
  2644
+        while True:
  2645
+            # apply all the tests
  2646
+            while l < base_len - 1 and\
  2647
+                  computed_words[l](base[l]) in orbit_reps[l] and\
  2648
+                  base_ordering[computed_words[l](base[l])] >\
  2649
+                  base_ordering[mu[l]] and\
  2650
+                  base_ordering[computed_words[l](base[l])] <\
  2651
+                  base_ordering[nu[l]] and\
  2652
+                  tests[l](computed_words[base_len - 1]):
  2653
+                # line 11: change the (partial) base of K
  2654
+                new_point = computed_words[l](base[l])
  2655
+                res_base[l] = new_point
  2656
+                temp_group = PermutationGroup(res_strong_gens_distr[l])
  2657
+                new_stab = temp_group.stabilizer(new_point)
  2658
+                res_strong_gens_distr[l + 1] = new_stab.generators
  2659
+                # line 12: calculate minimal orbit representatives for the
  2660
+                # l+1-th basic stabilizer
  2661
+                orbits = new_stab.orbits()
  2662
+                reps = []
  2663
+                for orbit in orbits:
  2664
+                    rep = min(orbit, key = lambda point: base_ordering[point])
  2665
+                    reps.append(rep)
  2666
+                orbit_reps[l + 1] = reps
  2667
+                # line 13: amend sorted orbits
  2668
+                l += 1
  2669
+                temp_orbit = [computed_words[l-1](point) for point\
  2670
+                             in basic_orbits[l]]
  2671
+                temp_orbit.sort(key = lambda point: base_ordering[point])
  2672
+                sorted_orbits[l] = temp_orbit
  2673
+                # lines 14 and 15: update variables used minimality tests
  2674
+                new_mu = degree + 1
  2675
+                for i in range(l):
  2676
+                    if base[l] in res_basic_orbits_init_base[i]:
  2677
+                        candidate = computed_words[i](base[i])
  2678
+                        if base_ordering[candidate] > base_ordering[new_mu]:
  2679
+                            new_mu = candidate
  2680
+                mu[l] = new_mu
  2681
+                temp_index = len(basic_orbits[l]) + 1 -\
  2682
+                             len(res_basic_orbits_init_base[l])
  2683
+                if temp_index >= len(sorted_orbits[l]):
  2684
+                    nu[l] = base_ordering[degree]
  2685
+                else:
  2686
+                    nu[l] = sorted_orbits[l][temp_index]
  2687
+                # line 16: determine the new transversal element
  2688
+                c[l] = 0
  2689
+                temp_point = sorted_orbits[l][c[l]]
  2690
+                temp_element = ~(computed_words[l - 1])
  2691
+                gamma = temp_element(temp_point)
  2692
+                u[l] = transversals[l][gamma]
  2693
+                # update computed words
  2694
+                computed_words[l] = computed_words[l-1] * u[l]
  2695
+            # lines 17 & 18: apply the tests to the group element found
  2696
+            g = computed_words[l]
  2697
+            temp_point = g(base[l])
  2698
+            if l == base_len - 1 and\
  2699
+               base_ordering[temp_point] > base_ordering[mu[l]] and\
  2700
+               base_ordering[temp_point] < base_ordering[nu[l]] and\
  2701
+               temp_point in orbit_reps[l] and\
  2702
+               tests[l](g) and\
  2703
+               prop(g):
  2704
+                # line 19: reset the base of K
  2705
+                gens = res.generators[:]
  2706
+                gens.append(g)
  2707
+                res = PermutationGroup(gens)
  2708
+                res_base = base[:]
  2709
+                # line 20: recalculate basic orbits (and transversals)
  2710
+                res_strong_gens.append(g)
  2711
+                res_strong_gens_distr = _distribute_gens_by_base(res_base,\
  2712
+                                                          res_strong_gens)
  2713
+                res_basic_orbits_init_base =\
  2714
+                [PermutationGroup(res_strong_gens_distr[i]).orbit(res_base[i])\
  2715
+                 for i in range(base_len)]
  2716
+                # line 21: recalculate orbit representatives
  2717
+                stab_f = PermutationGroup(res_strong_gens_distr[f])
  2718
+                temp_orbits = stab_f.orbits()
  2719
+                reps = []
  2720
+                for orbit in orbits:
  2721
+                    rep = min(orbit, key = lambda point: base_ordering[point])
  2722
+                    reps.append(rep)
  2723
+                orbit_reps[f] = reps
  2724
+                # line 22: reset the search depth
  2725
+                l = f
  2726
+            # line 23: go up the tree until in the first branch not fully
  2727
+            # searched
  2728
+            while l >= 0 and c[l] == len(basic_orbits[l]) - 1:
  2729
+                l = l - 1
  2730
+            # line 24: if the entire tree is traversed, return K
  2731
+            if l == -1:
  2732
+                return res
  2733
+            # lines 25-27: update orbit representatives
  2734
+            if l < f:
  2735
+                # line 26
  2736
+                f = l
  2737
+                c[l] = 0
  2738
+                # line 27
  2739
+                stab_f = PermutationGroup(res_strong_gens_distr[f])
  2740
+                temp_orbits = stab_f.orbits()
  2741
+                reps = []
  2742
+                for orbit in orbits:
  2743
+                    rep = min(orbit, key = lambda point: base_ordering[point])
  2744
+                    reps.append(rep)
  2745
+                orbit_reps[f] = reps
  2746
+                # line 28: update variables used for minimality testing
  2747
+                mu[l] = degree + 1
  2748
+                temp_index = len(basic_orbits[l]) + 1 -\
  2749
+                             len(res_basic_orbits_init_base[l])
  2750
+                if temp_index >= len(sorted_orbits[l]):
  2751
+                    nu[l] = base_ordering[degree]
  2752
+                else:
  2753
+                    nu[l] = sorted_orbits[l][temp_index]
  2754
+            # line 29: set the next element from the current branch and update
  2755
+            # accorndingly
  2756
+            c[l] += 1
  2757
+            element = ~(computed_words[l - 1])
  2758
+            gamma  = element(sorted_orbits[l][c[l]])
  2759
+            u[l] = transversals[l][gamma]
  2760
+            if l == 0:
  2761
+                computed_words[l] = u[l]
  2762
+            else:
  2763
+                computed_words[l] = computed_words[l - 1]*u[l]
  2764
+
2364 2765
     @property
2365 2766
     def transitivity_degree(self):
2366 2767
         """
55  sympy/combinatorics/tests/test_perm_groups.py
@@ -376,3 +376,58 @@ def test_direct_product_n():
376 376
     H = DirectProduct(D, C)
377 377
     assert H.order() == 32
378 378
     assert H.is_abelian == False
  379
+
  380
+def test_schreier_sims_incremental():
  381
+    identity = Permutation([0, 1, 2, 3, 4])
  382
+    TrivialGroup = PermutationGroup([identity])
  383
+    base, strong_gens = TrivialGroup.schreier_sims_incremental(base=[0, 1, 2])
  384
+    assert _verify_bsgs(TrivialGroup, base, strong_gens) == True
  385
+    S = SymmetricGroup(5)
  386
+    base, strong_gens = S.schreier_sims_incremental(base=[0,1,2])
  387
+    assert _verify_bsgs(S, base, strong_gens) == True
  388
+    D = DihedralGroup(2)
  389
+    base, strong_gens = D.schreier_sims_incremental(base=[1])
  390
+    assert _verify_bsgs(D, base, strong_gens) == True
  391
+    A = AlternatingGroup(7)
  392
+    gens = A.generators[:]
  393
+    gen0 = gens[0]
  394
+    gen1 = gens[1]
  395
+    gen1 = gen1*(~gen0)
  396
+    gen0 = gen0*gen1
  397
+    gen1 = gen0*gen1
  398
+    base, strong_gens = A.schreier_sims_incremental(base=[0,1], gens=gens)
  399
+    assert _verify_bsgs(A, base, strong_gens) == True
  400
+    C = CyclicGroup(11)
  401
+    gen = C.generators[0]
  402
+    base, strong_gens = C.schreier_sims_incremental(gens=[gen**3])
  403
+    assert _verify_bsgs(C, base, strong_gens) == True
  404
+
  405
+def test_subgroup_search():
  406
+    prop_true = lambda x: True
  407
+    prop_fix_points = lambda x: [x(point) for point in points] == points
  408
+    prop_comm_g = lambda x: x*g == g*x
  409
+    prop_even = lambda x: x.is_even
  410
+    for i in range(10, 17, 2):
  411
+        S = SymmetricGroup(i)
  412
+        A = AlternatingGroup(i)
  413
+        C = CyclicGroup(i)
  414
+        Sym = S.subgroup_search(prop_true)
  415
+        assert Sym == S
  416
+        Alt = S.subgroup_search(prop_even)
  417
+        assert Alt == A
  418
+        Sym = S.subgroup_search(prop_true, init_subgroup=C)
  419
+        assert Sym == S
  420
+        points = [7]
  421
+        assert S.stabilizer(7) == S.subgroup_search(prop_fix_points)
  422
+        points = [3, 4]
  423
+        assert S.stabilizer(3).stabilizer(4) == S.subgroup_search(prop_fix_points)
  424
+        points = [3, 5]
  425
+        fix35 = A.subgroup_search(prop_fix_points)
  426
+        points = [5]
  427
+        fix5 = A.subgroup_search(prop_fix_points)
  428
+        assert A.subgroup_search(prop_fix_points, init_subgroup=fix35) == fix5
  429
+        base, strong_gens = A.schreier_sims_incremental()
  430
+        g = A.generators[0]
  431
+        comm_g = A.subgroup_search(prop_comm_g, base=base, strong_gens=strong_gens)
  432
+        assert _verify_bsgs(comm_g, base, comm_g.generators) == True
  433
+        assert [prop_comm_g(gen) == True for gen in comm_g.generators]
28  sympy/combinatorics/tests/test_util.py
@@ -5,7 +5,7 @@
5 5
 from sympy.combinatorics.util import _check_cycles_alt_sym, _strip,\
6 6
 _distribute_gens_by_base, _strong_gens_from_distr,\
7 7
 _orbits_transversals_from_bsgs, _handle_precomputed_bsgs, _base_ordering,\
8  
-_verify_bsgs
  8
+_verify_bsgs, _remove_gens
9 9
 
10 10
 def test_check_cycles_alt_sym():
11 11
     perm1 = Permutation([[0, 1, 2, 3, 4, 5, 6], [7], [8], [9]])
@@ -44,9 +44,9 @@ def test_distribute_gens_by_base():
44 44
                                                    Permutation([0, 1, 3, 2])]]
45 45
 
46 46
 def test_strong_gens_from_distr():
47  
-    distr_gens = [[Permutation([0, 2, 1]), Permutation([1, 2, 0]),\
  47
+    strong_gens_distr = [[Permutation([0, 2, 1]), Permutation([1, 2, 0]),\
48 48
                   Permutation([1, 0, 2])], [Permutation([0, 2, 1])]]
49  
-    assert _strong_gens_from_distr(distr_gens) ==\
  49
+    assert _strong_gens_from_distr(strong_gens_distr) ==\
50 50
                                                      [Permutation([0, 2, 1]),\
51 51
                                                      Permutation([1, 2, 0]),\
52 52
                                                      Permutation([1, 0, 2])]
@@ -56,8 +56,8 @@ def test_orbits_transversals_from_bsgs():
56 56
     S.schreier_sims()
57 57
     base = S.base
58 58
     strong_gens = S.strong_gens
59  
-    distr_gens = _distribute_gens_by_base(base, strong_gens)
60  
-    result = _orbits_transversals_from_bsgs(base, distr_gens)
  59
+    strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
  60
+    result = _orbits_transversals_from_bsgs(base, strong_gens_distr)
61 61
     orbits = result[0]
62 62
     transversals = result[1]
63 63
     base_len = len(base)
@@ -77,8 +77,8 @@ def test_handle_precomputed_bsgs():
77 77
     base = A.base
78 78
     strong_gens = A.strong_gens
79 79
     result = _handle_precomputed_bsgs(base, strong_gens)
80  
-    distr_gens = _distribute_gens_by_base(base, strong_gens)
81  
-    assert distr_gens == result[2]
  80
+    strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
  81
+    assert strong_gens_distr == result[2]
82 82
     transversals = result[0]
83 83
     orbits = result[1]
84 84
     base_len = len(base)
@@ -106,3 +106,17 @@ def test_verify_bsgs():
106 106
     assert _verify_bsgs(S, base, strong_gens) == True
107 107
     assert _verify_bsgs(S, base[:-1], strong_gens) == False
108 108
     assert _verify_bsgs(S, base, S.generators) == False
  109
+
  110
+def test_remove_gens():
  111
+    S = SymmetricGroup(10)
  112
+    base, strong_gens = S.schreier_sims_incremental()
  113
+    new_gens = _remove_gens(base, strong_gens)
  114
+    assert _verify_bsgs(S, base, new_gens) == True
  115
+    A = AlternatingGroup(7)
  116
+    base, strong_gens = A.schreier_sims_incremental()
  117
+    new_gens = _remove_gens(base, strong_gens)
  118
+    assert _verify_bsgs(A, base, new_gens) == True
  119
+    D = DihedralGroup(2)
  120
+    base, strong_gens = D.schreier_sims_incremental()
  121
+    new_gens = _remove_gens(base, strong_gens)
  122
+    assert _verify_bsgs(D, base, new_gens) == True
125  sympy/combinatorics/util.py
@@ -172,7 +172,7 @@ def _distribute_gens_by_base(base, gens):
172 172
     return stabs
173 173
 
174 174
 def _handle_precomputed_bsgs(base, strong_gens, transversals=None,\
175  
-                             basic_orbits=None, distr_gens=None):
  175
+                             basic_orbits=None, strong_gens_distr=None):
176 176
     """
177 177
     Calculate BSGS-related structures from those present.
178 178
 
@@ -187,15 +187,15 @@ def _handle_precomputed_bsgs(base, strong_gens, transversals=None,\
187 187
     ``strong_gens`` - the strong generators
188 188
     ``transversals`` - basic transversals
189 189
     ``basic_orbits`` - basic orbits
190  
-    ``distr_gens`` - strong generators distributed by membership in basic
  190
+    ``strong_gens_distr`` - strong generators distributed by membership in basic
191 191
     stabilizers
192 192
 
193 193
     Returns
194 194
     =======
195 195
 
196  
-    ``(transversals, basic_orbits, distr_gens)`` where ``transversals`` are the
  196
+    ``(transversals, basic_orbits, strong_gens_distr)`` where ``transversals`` are the
197 197
     basic transversals, ``basic_orbits`` are the basic orbits, and
198  
-    ``distr_gens`` are the strong generators distributed by membership in basic
  198
+    ``strong_gens_distr`` are the strong generators distributed by membership in basic
199 199
     stabilizers.
200 200
 
201 201
     Examples
@@ -218,15 +218,15 @@ def _handle_precomputed_bsgs(base, strong_gens, transversals=None,\
218 218
     _orbits_transversals_from_bsgs, distribute_gens_by_base
219 219
 
220 220
     """
221  
-    if distr_gens is None:
222  
-        distr_gens = _distribute_gens_by_base(base, strong_gens)
  221
+    if strong_gens_distr is None:
  222
+        strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
223 223
     if transversals is None:
224 224
         if basic_orbits is None:
225 225
             basic_orbits, transversals =\
226  
-            _orbits_transversals_from_bsgs(base, distr_gens)
  226
+            _orbits_transversals_from_bsgs(base, strong_gens_distr)
227 227
         else:
228 228
             transversals =\
229  
-            _orbits_transversals_from_bsgs(base, distr_gens,
  229
+            _orbits_transversals_from_bsgs(base, strong_gens_distr,
230 230
                                            transversals_only=True)
231 231
     else:
232 232
         if basic_orbits is None:
@@ -234,9 +234,9 @@ def _handle_precomputed_bsgs(base, strong_gens, transversals=None,\
234 234
             basic_orbits = [None]*base_len
235 235
             for i in xrange(base_len):
236 236
                 basic_orbits[i] = transversals[i].keys()
237  
-    return transversals, basic_orbits, distr_gens
  237
+    return transversals, basic_orbits, strong_gens_distr
238 238
 
239  
-def _orbits_transversals_from_bsgs(base, distr_gens,\
  239
+def _orbits_transversals_from_bsgs(base, strong_gens_distr,\
240 240
                                    transversals_only=False):
241 241
     """
242 242
     Compute basic orbits and transversals from a base and strong generating set.
@@ -249,7 +249,7 @@ def _orbits_transversals_from_bsgs(base, distr_gens,\
249 249
     ==========
250 250
 
251 251
     ``base`` - the base
252  
-    ``distr_gens`` - strong generators distributed by membership in basic
  252
+    ``strong_gens_distr`` - strong generators distributed by membership in basic
253 253
     stabilizers
254 254
     ``transversals_only`` - a flag swithing between returning only the
255 255
     transversals/ both orbits and transversals
@@ -263,8 +263,8 @@ def _orbits_transversals_from_bsgs(base, distr_gens,\
263 263
     ... _distribute_gens_by_base)
264 264
     >>> S = SymmetricGroup(3)
265 265
     >>> S.schreier_sims()
266  
-    >>> distr_gens = _distribute_gens_by_base(S.base, S.strong_gens)
267  
-    >>> _orbits_transversals_from_bsgs(S.base, distr_gens)
  266
+    >>> strong_gens_distr = _distribute_gens_by_base(S.base, S.strong_gens)
  267
+    >>> _orbits_transversals_from_bsgs(S.base, strong_gens_distr)
268 268
     ([[0, 1, 2], [1, 2]], [{0: Permutation([0, 1, 2]),\
269 269
     1: Permutation([1, 2, 0]), 2: Permutation([2, 0, 1])},\
270 270
     {1: Permutation([0, 1, 2]), 2: Permutation([0, 2, 1])}])
@@ -281,7 +281,7 @@ def _orbits_transversals_from_bsgs(base, distr_gens,\
281 281
     if transversals_only is False:
282 282
         basic_orbits = [None]*base_len
283 283
     for i in xrange(base_len):
284  
-        group = PermutationGroup(distr_gens[i])
  284
+        group = PermutationGroup(strong_gens_distr[i])
285 285
         transversals[i] = dict(group.orbit_transversal(base[i], pairs=True))
286 286
         if transversals_only is False:
287 287
             basic_orbits[i] = transversals[i].keys()
@@ -290,6 +290,83 @@ def _orbits_transversals_from_bsgs(base, distr_gens,\
290 290
     else:
291 291
         return basic_orbits, transversals
292 292
 
  293
+def _remove_gens(base, strong_gens, basic_orbits=None, strong_gens_distr=None):
  294
+    """
  295
+    Remove redundant generators from a strong generating set.
  296
+
  297
+    Parameters
  298
+    ==========
  299
+
  300
+    ``base`` - a base
  301
+    ``strong_gens`` - a strong generating set relative to ``base``
  302
+    ``basic_orbits`` - basic orbits
  303
+    ``strong_gens_distr`` - strong generators distributed by membership in basic
  304
+    stabilizers
  305
+
  306
+    Returns
  307
+    =======
  308
+
  309
+    A strong generating set with respect to ``base`` which is a subset of
  310
+    ``strong_gens``.
  311
+
  312
+    Examples
  313
+    ========
  314
+
  315
+    >>> from sympy.combinatorics.named_groups import SymmetricGroup
  316
+    >>> from sympy.combinatorics.perm_groups import PermutationGroup
  317
+    >>> from sympy.combinatorics.util import _remove_gens, _verify_bsgs
  318
+    >>> S = SymmetricGroup(15)
  319
+    >>> base, strong_gens = S.schreier_sims_incremental()
  320
+    >>> len(strong_gens)
  321
+    26
  322
+    >>> new_gens = _remove_gens(base, strong_gens)
  323
+    >>> len(new_gens)
  324
+    14
  325
+    >>> _verify_bsgs(S, base, new_gens)
  326
+    True
  327
+
  328
+    Notes
  329
+    =====
  330
+
  331
+    This procedure is outlined in [1],p.95.
  332
+
  333
+    References
  334
+    ==========
  335
+
  336
+    [1] Holt, D., Eick, B., O'Brien, E.
  337
+    "Handbook of computational group theory"
  338
+
  339
+    """
  340
+    from sympy.combinatorics.perm_groups import PermutationGroup
  341
+    base_len = len(base)
  342
+    degree = strong_gens[0].size
  343
+    identity = _new_from_array_form(range(degree))
  344
+    if strong_gens_distr is None:
  345
+        strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
  346
+    temp = strong_gens_distr[:]
  347
+    if basic_orbits is None:
  348
+        basic_orbits = []
  349
+        for i in range(base_len):
  350
+            stab = PermutationGroup(strong_gens_distr[i])
  351
+            basic_orbit = stab.orbit(base[i])
  352
+            basic_orbits.append(basic_orbit)
  353
+    strong_gens_distr.append([])
  354
+    res = strong_gens[:]
  355
+    for i in range(base_len - 1, -1, -1):
  356
+        gens_copy = strong_gens_distr[i][:]
  357
+        for gen in strong_gens_distr[i]:
  358
+            if gen not in strong_gens_distr[i + 1]:
  359
+                temp_gens = gens_copy[:]
  360
+                temp_gens.remove(gen)
  361
+                if temp_gens == []:
  362
+                    continue
  363
+                temp_group = PermutationGroup(temp_gens)
  364
+                temp_orbit = temp_group.orbit(base[i])
  365
+                if temp_orbit == basic_orbits[i]:
  366
+                    gens_copy.remove(gen)
  367
+                    res.remove(gen)
  368
+    return res
  369
+
293 370
 def _strip(g, base, orbs, transversals):
294 371
     """
295 372
     Attempt to decompose a permutation using a (possibly partial) BSGS
@@ -363,7 +440,7 @@ def _strip(g, base, orbs, transversals):
363 440
         h = ~u*h
364 441
     return h, base_len + 1
365 442
 
366  
-def _strong_gens_from_distr(distr_gens):
  443
+def _strong_gens_from_distr(strong_gens_distr):
367 444
     """
368 445
     Retrieve strong generating set from generators of basic stabilizers.
369 446
 
@@ -373,7 +450,7 @@ def _strong_gens_from_distr(distr_gens):
373 450
     Parameters
374 451
     ==========
375 452
 
376  
-    ``distr_gens`` - strong generators distributed by membership in basic
  453
+    ``strong_gens_distr`` - strong generators distributed by membership in basic
377 454
     stabilizers
378 455
 
379 456
     Examples
@@ -386,8 +463,8 @@ def _strong_gens_from_distr(distr_gens):
386 463
     >>> S.schreier_sims()
387 464
     >>> S.strong_gens
388 465
     [Permutation([1, 2, 0]), Permutation([1, 0, 2]), Permutation([0, 2, 1])]
389  
-    >>> distr_gens = _distribute_gens_by_base(S.base, S.strong_gens)
390  
-    >>> _strong_gens_from_distr(distr_gens)
  466
+    >>> strong_gens_distr = _distribute_gens_by_base(S.base, S.strong_gens)
  467
+    >>> _strong_gens_from_distr(strong_gens_distr)
391 468
     [Permutation([1, 2, 0]), Permutation([1, 0, 2]), Permutation([0, 2, 1])]
392 469
 
393 470