diff --git a/doc/src/modules/combinatorics/group_constructs.rst b/doc/src/modules/combinatorics/group_constructs.rst
new file mode 100644
index 000000000000..3861622cdd3c
--- /dev/null
+++ b/doc/src/modules/combinatorics/group_constructs.rst
@@ -0,0 +1,8 @@
+.. _combinatorics-group_constructs:
+
+Group constructors
+==================
+
+.. module:: sympy.combinatorics.group_constructs
+
+.. autofunction:: DirectProduct
diff --git a/doc/src/modules/combinatorics/index.rst b/doc/src/modules/combinatorics/index.rst
index 9891816d1ab0..617195b0cdaa 100644
--- a/doc/src/modules/combinatorics/index.rst
+++ b/doc/src/modules/combinatorics/index.rst
@@ -17,3 +17,5 @@ Contents
    graycode.rst
    named_groups.rst
    util.rst
+   group_constructs.rst
+   testutil.rst
diff --git a/doc/src/modules/combinatorics/testutil.rst b/doc/src/modules/combinatorics/testutil.rst
new file mode 100644
index 000000000000..507510b9a2de
--- /dev/null
+++ b/doc/src/modules/combinatorics/testutil.rst
@@ -0,0 +1,16 @@
+.. _combinatorics-testutil:
+
+Test Utilities
+==============
+
+.. module:: sympy.combinatorics.testutil
+
+.. autofunction:: _cmp_perm_lists
+
+.. autofunction:: _naive_list_centralizer
+
+.. autofunction:: _verify_bsgs
+
+.. autofunction:: _verify_centralizer
+
+.. autofunction:: _verify_normal_closure
diff --git a/doc/src/modules/combinatorics/util.rst b/doc/src/modules/combinatorics/util.rst
index 2901a8e4505a..2c49bcbac10c 100644
--- a/doc/src/modules/combinatorics/util.rst
+++ b/doc/src/modules/combinatorics/util.rst
@@ -20,5 +20,3 @@ Utilities
 .. autofunction:: _strip
 
 .. autofunction:: _strong_gens_from_distr
-
-.. autofunction:: _verify_bsgs
diff --git a/sympy/combinatorics/group_constructs.py b/sympy/combinatorics/group_constructs.py
new file mode 100644
index 000000000000..0dcbd41cb771
--- /dev/null
+++ b/sympy/combinatorics/group_constructs.py
@@ -0,0 +1,53 @@
+from sympy.combinatorics.perm_groups import PermutationGroup
+from sympy.combinatorics.permutations import _new_from_array_form
+
+def DirectProduct(*groups):
+    """
+    Returns the direct product of several groups as a permutation group.
+
+    This is implemented much like the __mul__ procedure for taking the direct
+    product of two permutation groups, but the idea of shifting the
+    generators is realized in the case of an arbitrary number of groups.
+    A call to DirectProduct(G1, G2, ..., Gn) is generally expected to be faster
+    than a call to G1*G2*...*Gn (and thus the need for this algorithm).
+
+    Examples
+    ========
+
+    >>> from sympy.combinatorics.group_constructs import DirectProduct
+    >>> from sympy.combinatorics.named_groups import CyclicGroup
+    >>> C = CyclicGroup(4)
+    >>> G = DirectProduct(C,C,C)
+    >>> G.order()
+    64
+
+    See Also
+    ========
+    __mul__
+
+    """
+    degrees = []
+    gens_count = []
+    total_degree = 0
+    total_gens = 0
+    for group in groups:
+        current_deg = group.degree
+        current_num_gens = len(group.generators)
+        degrees.append(current_deg)
+        total_degree += current_deg
+        gens_count.append(current_num_gens)
+        total_gens += current_num_gens
+    array_gens = []
+    for i in range(total_gens):
+        array_gens.append(range(total_degree))
+    current_gen = 0
+    current_deg = 0
+    for i in xrange(len(gens_count)):
+        for j in xrange(current_gen, current_gen + gens_count[i]):
+            gen = ((groups[i].generators)[j - current_gen]).array_form
+            array_gens[j][current_deg:current_deg + degrees[i]] =\
+            [ x + current_deg for x in gen]
+        current_gen += gens_count[i]
+        current_deg += degrees[i]
+    perm_gens = [_new_from_array_form(array) for array in array_gens]
+    return PermutationGroup(perm_gens)
diff --git a/sympy/combinatorics/named_groups.py b/sympy/combinatorics/named_groups.py
index 251366946897..99f86a31a896 100644
--- a/sympy/combinatorics/named_groups.py
+++ b/sympy/combinatorics/named_groups.py
@@ -1,4 +1,5 @@
-from sympy.combinatorics.perm_groups import PermutationGroup, DirectProduct
+from sympy.combinatorics.perm_groups import PermutationGroup
+from sympy.combinatorics.group_constructs import DirectProduct
 from sympy.combinatorics.permutations import Permutation, _new_from_array_form
 
 def AbelianGroup(*cyclic_orders):
diff --git a/sympy/combinatorics/perm_groups.py b/sympy/combinatorics/perm_groups.py
index d476879a6d07..2315cdd203f9 100644
--- a/sympy/combinatorics/perm_groups.py
+++ b/sympy/combinatorics/perm_groups.py
@@ -340,6 +340,14 @@ class PermutationGroup(Basic):
 
     [8] http://en.wikipedia.org/wiki/Multiply_transitive_group#Multiply_transitive_groups
 
+    [9] http://en.wikipedia.org/wiki/Center_%28group_theory%29
+
+    [10] http://en.wikipedia.org/wiki/Centralizer_and_normalizer
+
+    [11] http://groupprops.subwiki.org/wiki/Derived_subgroup
+
+    [12] http://en.wikipedia.org/wiki/Nilpotent_group
+
     """
 
     def __eq__(self, gr):
@@ -423,6 +431,9 @@ def __new__(cls, *args, **kw_args):
         obj._is_sym = None
         obj._is_alt = None
         obj._is_primitive = None
+        obj._is_nilpotent = None
+        obj._is_solvable = None
+        obj._is_trivial = None
         obj._transitivity_degree = None
         obj._max_div = None
         size = len(args[0][0].array_form)
@@ -639,8 +650,8 @@ def baseswap(self, base, strong_gens, pos, randomized=False,\
         ``randomized`` - switch between randomized and deterministic version
         ``transversals`` - transversals for the basic orbits, if known
         ``basic_orbits`` - basic orbits, if known
-        ``strong_gens_distr`` - strong generators distributed by basic stabilizers,
-        if known
+        ``strong_gens_distr`` - strong generators distributed by basic
+        stabilizers, if known
 
         Returns
         =======
@@ -652,7 +663,7 @@ def baseswap(self, base, strong_gens, pos, randomized=False,\
         ========
 
         >>> from sympy.combinatorics.named_groups import SymmetricGroup
-        >>> from sympy.combinatorics.util import _verify_bsgs
+        >>> from sympy.combinatorics.testutil import _verify_bsgs
         >>> S = SymmetricGroup(4)
         >>> S.schreier_sims()
         >>> S.baseswap(S.base, S.strong_gens, 1, randomized=False)
@@ -742,6 +753,7 @@ def baseswap(self, base, strong_gens, pos, randomized=False,\
             if gen not in strong_gens_new:
                 strong_gens_new.append(gen)
         return base_new, strong_gens_new
+
     @property
     def basic_orbits(self):
         """
@@ -841,46 +853,194 @@ def basic_transversals(self):
             self.schreier_sims()
         return self._transversals
 
-    def commutator(self):
+    def center(self):
+        r"""
+        Return the center of a permutation group.
+
+        The center for a group `G` is defined as
+        `Z(G) = \{z\in G | \forall g\in G, zg = gz \}`,
+        the set of elements of `G` that commute with all elements of `G`.
+        It is equal to the centralizer of `G` inside `G`, and is naturally a
+        subgroup of `G` ([9]).
+
+        Examples
+        ========
+
+        >>> from sympy.combinatorics.perm_groups import PermutationGroup
+        >>> from sympy.combinatorics.named_groups import DihedralGroup
+        >>> D = DihedralGroup(4)
+        >>> G = D.center()
+        >>> G.order()
+        2
+
+        See Also
+        ========
+
+        centralizer
+
+        Notes
+        =====
+
+        This is a naive implementation that is a straightforward application
+        of ``.centralizer()``
+
         """
-        commutator subgroup
+        return self.centralizer(self)
 
-        The commutator subgroup is the subgroup generated by all
-        commutators; it is equal to the normal closure of the set
-        of commutators of the generators.
+    def centralizer(self, other):
+        r"""
+        Return the centralizer of a group/set/element.
+
+        The centralizer of a set of permutations `S` inside
+        a group `G` is the set of elements of `G` that commute with all
+        elements of `S`:
+        `C_G(S) = \{ g \in G | gs = sg \forall s \in S\}` ([10])
+        Usually, `S` is a subset of `G`, but if `G` is a proper subgroup of
+        the full symmetric group, we allow for `S` to have elements outside
+        `G`.
+        It is naturally a subgroup of `G`; the centralizer of a permutation
+        group is equal to the centralizer of any set of generators for that
+        group, since any element commuting with the generators commutes with
+        any product of the  generators.
 
-        see http://groupprops.subwiki.org/wiki/Derived_subgroup
+        Parameters
+        ==========
+
+        ``other`` - a permutation group/list of permutations/single permutation
 
         Examples
         ========
 
-        >>> from sympy.combinatorics.permutations import Permutation
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
-        >>> a = Permutation([1, 0, 2, 4, 3])
-        >>> b = Permutation([0, 1, 3, 2, 4])
-        >>> G = PermutationGroup([a, b])
-        >>> C = G.commutator()
-        >>> list(C.generate(af=True))
-        [[0, 1, 2, 3, 4], [0, 1, 3, 4, 2], [0, 1, 4, 2, 3]]
+        >>> from sympy.combinatorics.named_groups import (SymmetricGroup,
+        ... CyclicGroup)
+        >>> S = SymmetricGroup(6)
+        >>> C = CyclicGroup(6)
+        >>> H = S.centralizer(C)
+        >>> H == C
+        True
+
+        See Also
+        ========
+
+        subgroup_search
+
+        Notes
+        =====
+
+        The implementation is an application of ``.subgroup_search()`` with
+        tests using a specific base for the group `G`.
+
+        """
+        if hasattr(other, 'generators'):
+            if other.is_trivial or self.is_trivial:
+                return self
+            degree = self.degree
+            identity = _new_from_array_form(range(degree))
+            orbits = other.orbits()
+            num_orbits = len(orbits)
+            orbits.sort(key = lambda x: -len(x))
+            long_base = []
+            orbit_reps = [None]*num_orbits
+            orbit_reps_indices = [None]*num_orbits
+            orbit_descr = [None]*degree
+            for i in range(num_orbits):
+                orbit = list(orbits[i])
+                orbit_reps[i] = orbit[0]
+                orbit_reps_indices[i] = len(long_base)
+                for point in orbit:
+                    orbit_descr[point] = i
+                long_base = long_base + orbit
+            base, strong_gens = self.schreier_sims_incremental(base=long_base)
+            strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
+            i = 0
+            for i in range(len(base)):
+                if strong_gens_distr[i] == [identity]:
+                    break
+            base = base[:i]
+            base_len = i
+            for j in range(num_orbits):
+                if base[base_len - 1] in orbits[j]:
+                    break
+            rel_orbits = orbits[: j + 1]
+            num_rel_orbits = len(rel_orbits)
+            transversals = [None]*num_rel_orbits
+            for j in range(num_rel_orbits):
+                rep = orbit_reps[j]
+                transversals[j] = dict(
+                                      other.orbit_transversal(rep, pairs=True))
+            trivial_test = lambda x: True
+            tests = [None]*base_len
+            for l in range(base_len):
+                if base[l] in orbit_reps:
+                    tests[l] = trivial_test
+                else:
+                    def test(computed_words, l=l):
+                        g = computed_words[l]
+                        rep_orb_index = orbit_descr[base[l]]
+                        rep = orbit_reps[rep_orb_index]
+                        rep_orb = orbits[rep_orb_index]
+                        im = g(base[l])
+                        im_rep = g(rep)
+                        tr_el = transversals[rep_orb_index][base[l]]
+                        if im != tr_el(im_rep):
+                            return False
+                        else:
+                            return True
+                    tests[l] = test
+            def prop(g):
+                return [g*gen for gen in other.generators] ==\
+                       [gen*g for gen in other.generators]
+            return self.subgroup_search(prop, base=base,
+                                        strong_gens=strong_gens, tests=tests)
+        elif hasattr(other, '__getitem__'):
+            gens = list(other)
+            return self.centralizer(PermutationGroup(gens))
+        elif hasattr(other, 'array_form'):
+            return self.centralizer(PermutationGroup([other]))
+
+    def commutator(self, G, H):
+        """
+        Return the commutator of two subgroups.
+
+        For a permutation group `K` and subgroups `G`, `H`, the
+        commutator of `G` and `H` is defined as the group generated
+        by all the commutators `[g, h] = hgh^{-1}g^{-1}` for `g` in `G` and
+        `h` in `H`. It is naturally a subgroup of `K` ([1],p.27).
+
+        Examples
+        ========
+
+        >>> from sympy.combinatorics.named_groups import (SymmetricGroup,
+        ... AlternatingGroup)
+        >>> S = SymmetricGroup(5)
+        >>> A = AlternatingGroup(5)
+        >>> G = S.commutator(S, A)
+        >>> G == A
+        True
+
+        See Also
+        ========
+
+        derived_subgroup
+
+        Notes
+        =====
+
+        The commutator of two subgroups `H, G` is equal to the normal closure
+        of the commutators of all the generators, i.e. `hgh^{-1}g^{-1}` for `h`
+        a generator of `H` and `g` a generator of `G` ([1],p.28)
 
         """
-        r = self._r
-        gens = [p.array_form for p in self.generators]
-        gens_inv = [perm_af_invert(p) for p in gens]
-        set_commutators = set()
-        for i in range(r):
-            for j in range(r):
-                p1 = gens[i]
-                p1inv = gens_inv[i]
-                p2 = gens[j]
-                p2inv = gens_inv[j]
-                c = [p1[p2[p1inv[k]]] for k in p2inv]
-                ct = tuple(c)
-                if not ct in set_commutators:
-                    set_commutators.add(ct)
-        cms = [Permutation(p) for p in set_commutators]
-        G2 = self.normal_closure(cms)
-        return G2
+        ggens = G.generators
+        hgens = H.generators
+        commutators = []
+        for ggen in ggens:
+            for hgen in hgens:
+                commutator = hgen*ggen*(~hgen)*(~ggen)
+                if commutator not in commutators:
+                    commutators.append(commutator)
+        res = self.normal_closure(commutators)
+        return res
 
     def coset_decomposition(self, g):
         """
@@ -1061,6 +1221,98 @@ def degree(self):
         """
         return self._degree
 
+    def derived_series(self):
+        r"""
+        Return the derived series for the group.
+
+        The derived series for a group `G` is defined as
+        `G = G_0 > G_1 > G_2 > \ldots`
+        where `G_i = [G_{i-1}, G_{i-1}]`, i.e. `G_i` is the derived subgroup of
+        `G_{i-1}`, for `i\in\mathbb{N}`. When we have `G_k = G_{k-1}` for some
+        `k\in\mathbb{N}`, the series terminates.
+
+        Returns
+        =======
+
+        A list of permutation groups containing the members of the derived
+        series in the order `G = G_0, G_1, G_2, \ldots`.
+
+        Examples
+        ========
+
+        >>> from sympy.combinatorics.named_groups import (SymmetricGroup,
+        ... AlternatingGroup, DihedralGroup)
+        >>> A = AlternatingGroup(5)
+        >>> len(A.derived_series())
+        1
+        >>> S = SymmetricGroup(4)
+        >>> len(S.derived_series())
+        4
+        >>> S.derived_series()[1] == AlternatingGroup(4)
+        True
+        >>> S.derived_series()[2] == DihedralGroup(2)
+        True
+
+        See Also
+        ========
+
+        derived_subgroup
+
+        """
+        res = [self]
+        current = self
+        next = self.derived_subgroup()
+        while current != next:
+            res.append(next)
+            current = next
+            next = next.derived_subgroup()
+        return res
+
+    def derived_subgroup(self):
+        """
+        Compute the derived subgroup.
+
+        The derived subgroup, or commutator subgroup is the subgroup generated
+        by all commutators `[g, h] = hgh^{-1}g^{-1}` for `g, h\in G` ; it is
+        equal to the normal closure of the set of commutators of the generators
+        ([1],p.28, [11]).
+
+        Examples
+        ========
+
+        >>> from sympy.combinatorics.permutations import Permutation
+        >>> from sympy.combinatorics.perm_groups import PermutationGroup
+        >>> a = Permutation([1, 0, 2, 4, 3])
+        >>> b = Permutation([0, 1, 3, 2, 4])
+        >>> G = PermutationGroup([a, b])
+        >>> C = G.derived_subgroup()
+        >>> list(C.generate(af=True))
+        [[0, 1, 2, 3, 4], [0, 1, 3, 4, 2], [0, 1, 4, 2, 3]]
+
+        See Also
+        ========
+
+        derived_series
+
+        """
+        r = self._r
+        gens = [p.array_form for p in self.generators]
+        gens_inv = [perm_af_invert(p) for p in gens]
+        set_commutators = set()
+        for i in range(r):
+            for j in range(r):
+                p1 = gens[i]
+                p1inv = gens_inv[i]
+                p2 = gens[j]
+                p2inv = gens_inv[j]
+                c = [p1[p2[p1inv[k]]] for k in p2inv]
+                ct = tuple(c)
+                if not ct in set_commutators:
+                    set_commutators.add(ct)
+        cms = [Permutation(p) for p in set_commutators]
+        G2 = self.normal_closure(cms)
+        return G2
+
     def generate(self, method="coset", af=False):
         """
         return iterator to generate the elements of the group
@@ -1281,8 +1533,8 @@ def is_alt_sym(self, eps=0.05, _random_prec=None):
         Monte Carlo test for the symmetric/alternating group for degrees >= 8.
 
         More specifically, it is one-sided Monte Carlo with the
-        answer True (i.e., G is symmetric/alternating) guaranteed to be correct,
-        and the answer False being incorrect with probability eps.
+        answer True (i.e., G is symmetric/alternating) guaranteed to be
+        correct, and the answer False being incorrect with probability eps.
 
         Notes
         =====
@@ -1337,6 +1589,50 @@ def is_alt_sym(self, eps=0.05, _random_prec=None):
                     return True
             return False
 
+    @property
+    def is_nilpotent(self):
+        """
+        Test if the group is nilpotent.
+
+        A group `G` is nilpotent if it has a central series of finite length.
+        Alternatively, `G` is nilpotent if its lower central series terminates
+        with the trivial group. Every nilpotent group is also solvable
+        ([1],p.29, [12]).
+
+        Examples
+        ========
+
+        >>> from sympy.combinatorics.named_groups import (SymmetricGroup,
+        ... CyclicGroup)
+        >>> C = CyclicGroup(6)
+        >>> C.is_nilpotent
+        True
+        >>> S = SymmetricGroup(5)
+        >>> S.is_nilpotent
+        False
+
+        See Also
+        ========
+
+        lower_central_series, is_solvable
+
+        """
+        if self._is_nilpotent is None:
+            lcs = self.lower_central_series()
+            terminator = lcs[len(lcs)-1]
+            gens = terminator.generators
+            degree = self.degree
+            identity = _new_from_array_form(range(degree))
+            if [identity for gen in gens] == gens:
+                self._is_solvable = True
+                self._is_nilpotent = True
+                return True
+            else:
+                self._is_nilpotent = False
+                return False
+        else:
+            return self._is_nilpotent
+
     def is_normal(self, gr):
         """
         test if G=self is a normal subgroup of gr
@@ -1428,38 +1724,42 @@ def is_primitive(self, randomized=True):
         self._is_primitive = True
         return True
 
+    @property
     def is_solvable(self):
         """
-        test if the group G is solvable
+        Test if the group  is solvable
 
-        G is solvable if the derived series
-        G = G_0 < G_1 < ... < G_k = 1, with G_{i+1} = G.commutator()
-        see http://en.wikipedia.org/wiki/Solvable_group
+        `G` is solvable if its derived series terminates with the trivial
+        group ([1],p.29).
 
         Examples
         ========
 
-        >>> from sympy.combinatorics.permutations import Permutation
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
-        >>> a = Permutation([1,2,0])
-        >>> b = Permutation([1,0,2])
-        >>> G = PermutationGroup([a, b])
-        >>> G.is_solvable()
+        >>> from sympy.combinatorics.named_groups import SymmetricGroup
+        >>> S = SymmetricGroup(3)
+        >>> S.is_solvable
         True
 
-        """
-        order = self.order()
-        if order == 1:
-            return True
-        G = self
-        while order > 1:
+        See Also
+        ========
+
+        is_nilpotent, derived_series
 
-            G = G.commutator()
-            order1 = G.order()
-            if order1 == order:
+        """
+        if self._is_solvable is None:
+            ds = self.derived_series()
+            terminator = ds[len(ds) - 1]
+            gens = terminator.generators
+            degree = self.degree
+            identity = _new_from_array_form(range(degree))
+            if [identity for gen in gens] == gens:
+                self._is_solvable = True
+                return True
+            else:
+                self._is_solvable = False
                 return False
-            order = order1
-        return True
+        else:
+            return self._is_solvable
 
     def is_subgroup(self, gr):
         """
@@ -1519,6 +1819,75 @@ def is_transitive(self):
         self._is_transitive = ans
         return ans
 
+    @property
+    def is_trivial(self):
+        """
+        Test if the group is the trivial group.
+
+        This is true if and only if all the generators are the identity.
+
+        Examples
+        ========
+
+        >>> from sympy.combinatorics.perm_groups import PermutationGroup
+        >>> from sympy.combinatorics.permutations import Permutation
+        >>> id = Permutation(range(5))
+        >>> G = PermutationGroup([id, id, id])
+        >>> G.is_trivial
+        True
+
+        """
+        if self._is_trivial is None:
+            gens = self.generators
+            degree = self.degree
+            identity = _new_from_array_form(range(degree))
+            res = [identity for gen in gens] == gens
+            self._is_trivial = res
+            return res
+        else:
+            return self._is_trivial
+
+    def lower_central_series(self):
+        r"""
+        Return the lower central series for the group.
+
+        The lower central series for a group `G` is the series
+        `G = G_0 > G_1 > G_2 > \ldots` where
+        `G_k = [G, G_{k-1}]`, i.e. every term after the first is equal to the
+        commutator of `G` and the previous term in `G1` ([1],p.29).
+
+        Returns
+        =======
+
+        A list of permutation groups in the order
+        `G = G_0, G_1, G_2, \ldots`
+
+        Examples
+        ========
+
+        >>> from sympy.combinatorics.named_groups import (AlternatingGroup,
+        ... DihedralGroup)
+        >>> A = AlternatingGroup(4)
+        >>> len(A.lower_central_series())
+        2
+        >>> A.lower_central_series()[1] == DihedralGroup(2)
+        True
+
+        See Also
+        ========
+
+        commutator, derived_series
+
+        """
+        res = [self]
+        current = self
+        next = self.commutator(self, current)
+        while current != next:
+            res.append(next)
+            current = next
+            next = self.commutator(self, current)
+        return res
+
     @property
     def max_div(self):
         """
@@ -1640,46 +2009,103 @@ def minimal_block(self, points):
             self._union_find_rep(i, parents)
         return parents
 
-    def normal_closure(self, gens):
-        """
-        normal closure in self of a list gens2 of permutations
+    def normal_closure(self, other, k=10):
+        r"""
+        Return the normal closure of a subgroup/set of permutations.
+
+        If `S` is a subset of a group `G`, the normal closure of `A` in `G`
+        is defined as the intersection of all normal subgroups of `G` that
+        contain `A` ([1],p.14). Alternatively, it is the group generated by
+        the conjugates `x^{-1}yx` for `x` a generator of `G` and `y` a
+        generator of the subgroup `\left\langle S\right\rangle` generated by
+        `S` (for some chosen generating set for `\left\langle S\right\rangle`)
+        ([1],p.73).
+
+        Parameters
+        ==========
+
+        ``other`` - a subgroup/list of permutations/single permutation
+        ``k`` - an implementation-specific parameter that determines the number
+        of conjugates that are adjoined to ``other`` at once
 
         Examples
         ========
 
-        >>> from sympy.combinatorics.permutations import Permutation
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
-        >>> a = Permutation([1, 2, 0])
-        >>> b = Permutation([1, 0, 2])
-        >>> G = PermutationGroup([a, b])
+        >>> from sympy.combinatorics.named_groups import (SymmetricGroup,
+        ... CyclicGroup, AlternatingGroup)
+        >>> S = SymmetricGroup(5)
+        >>> C = CyclicGroup(5)
+        >>> G = S.normal_closure(C)
         >>> G.order()
-        6
-        >>> G1 = G.normal_closure([a])
-        >>> list(G1.generate(af=True))
-        [[0, 1, 2], [1, 2, 0], [2, 0, 1]]
-
-        """
-        G2 = PermutationGroup(gens)
-        if G2.is_normal(self):
-            return G2
-        gens1 = [p.array_form for p in self.generators]
-        b = 0
-        while not b:
-            gens2 = [p.array_form for p in G2.generators]
-            b = 1
-            for g1 in gens1:
-                if not b:
-                    break
-                for g2 in gens2:
-                    p = perm_af_muln(g1, g2, perm_af_invert(g1))
-                    p = Permutation(p)
-                    if not G2.has_element(p):
-                        gens2 = G2.generators + [p]
-                        G2 = PermutationGroup(gens2)
-                        b = 0
-                        break
+        60
+        >>> G == AlternatingGroup(5)
+        True
 
-        return G2
+        See Also
+        ========
+
+        commutator, derived_subgroup, random_pr
+
+        Notes
+        =====
+
+        The algorithm is described in [1],pp.73-74; it makes use of the
+        generation of random elements for permutation groups by the product
+        replacement algorithm.
+
+        """
+        if hasattr(other, 'generators'):
+            degree = self.degree
+            identity = _new_from_array_form(range(degree))
+
+            if other.generators == [identity for gen in other.generators]:
+                return other
+            Z = PermutationGroup(other.generators[:])
+            base, strong_gens = Z.schreier_sims_incremental()
+            strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
+            basic_orbits, basic_transversals =\
+                    _orbits_transversals_from_bsgs(base, strong_gens_distr)
+            C = False
+            self._random_pr_init(r=10, n=20)
+
+            while C == False:
+                Z._random_pr_init(r=10, n=10)
+                for i in range(k):
+                    g = self.random_pr()
+                    h = Z.random_pr()
+                    conj = (~g)*h*g
+                    res = _strip(conj, base, basic_orbits, basic_transversals)
+                    if res[0] != identity or res[1] != len(base) + 1:
+                        gens = Z.generators
+                        gens.append(conj)
+                        Z = PermutationGroup(gens)
+                        strong_gens.append(conj)
+                        temp_base, temp_strong_gens =\
+                                Z.schreier_sims_incremental(base, strong_gens)
+                        base, strong_gens = temp_base, temp_strong_gens
+                        strong_gens_distr =\
+                                _distribute_gens_by_base(base, strong_gens)
+                        basic_orbits, basic_transversals =\
+                                _orbits_transversals_from_bsgs(base,\
+                                strong_gens_distr)
+                C = True
+                break_flag = False
+                for g in self.generators:
+                    for h in Z.generators:
+                        conj = (~g)*h*g
+                        res = _strip(conj, base, basic_orbits,\
+                                     basic_transversals)
+                        if res[0] != identity or res[1] != len(base) + 1:
+                            C = False
+                            break_flag = True
+                            break
+                    if break_flag == True:
+                        break
+            return Z
+        elif hasattr(other, '__getitem__'):
+            return self.normal_closure(PermutationGroup(other))
+        elif hasattr(other, 'array_form'):
+            return self.normal_closure(PermutationGroup([other]))
 
     def orbit(self, alpha, action='tuples'):
         r"""
@@ -1902,6 +2328,48 @@ def order(self):
             m *= x
         return m
 
+    def pointwise_stabilizer(self, points):
+        r"""
+        Return the pointwise stabilizer for a set of points.
+
+        For a permutation group `G` and a set of points
+        `\{p_1, p_2,\ldots, p_k\}`, the pointwise stabilizer of
+        `p_1, p_2, \ldots, p_k` is defined as
+        `G_{p_1,\ldots, p_k} =
+        \{g\in G | g(p_i) = p_i \forall i\in\{1, 2,\ldots,k\}\} ([1],p20).
+        It is a subgroup of `G`.
+
+        Examples
+        ========
+
+        >>> from sympy.combinatorics.named_groups import SymmetricGroup
+        >>> S = SymmetricGroup(7)
+        >>> Stab = S.pointwise_stabilizer([2, 3, 5])
+        >>> Stab == S.stabilizer(2).stabilizer(3).stabilizer(5)
+        True
+
+        See Also
+        ========
+
+        stabilizer, schreier_sims_incremental
+
+        Notes
+        =====
+
+        Rather than the obvious implementation using successive calls to
+        .stabilizer(), this uses the incremental Schreier-Sims algorithm
+        to obtain a base with starting segment - the given points.
+
+        """
+        base, strong_gens = self.schreier_sims_incremental(base=points)
+        stab_gens = []
+        degree = self.degree
+        identity = _new_from_array_form(range(degree))
+        for gen in strong_gens:
+            if [gen(point) for point in points] == points:
+                stab_gens.append(gen)
+        return PermutationGroup(stab_gens)
+
     def random(self, af=False):
         """
         return a random group element
@@ -2107,7 +2575,7 @@ def schreier_sims_incremental(self, base=None, gens=None):
 
         >>> from sympy.combinatorics.named_groups import AlternatingGroup
         >>> from sympy.combinatorics.perm_groups import PermutationGroup
-        >>> from sympy.combinatorics.util import _verify_bsgs
+        >>> from sympy.combinatorics.testutil import _verify_bsgs
         >>> A = AlternatingGroup(7)
         >>> base = [2, 3]
         >>> seq = [2, 3]
@@ -2202,7 +2670,8 @@ def schreier_sims_incremental(self, base=None, gens=None):
                             # data structures and start over
                             for l in range(i + 1, j):
                                 strong_gens_distr[l].append(h)
-                                stabs[l] = PermutationGroup(strong_gens_distr[l])
+                                stabs[l] =\
+                                        PermutationGroup(strong_gens_distr[l])
                                 transversals[l] =\
                                 dict(stabs[l].orbit_transversal(_base[l],\
                                                                 pairs=True))
@@ -2253,7 +2722,7 @@ def schreier_sims_random(self, base=None, gens=None, consec_succ=10,\
         ========
 
         >>> from sympy.combinatorics.perm_groups import PermutationGroup
-        >>> from sympy.combinatorics.util import _verify_bsgs
+        >>> from sympy.combinatorics.testutil import _verify_bsgs
         >>> from sympy.combinatorics.named_groups import SymmetricGroup
         >>> S = SymmetricGroup(5)
         >>> base, strong_gens = S.schreier_sims_random(consec_succ=5)
@@ -2541,7 +3010,7 @@ def subgroup_search(self, prop, base=None, strong_gens=None, tests=None,\
         >>> from sympy.combinatorics.named_groups import (SymmetricGroup,
         ... AlternatingGroup)
         >>> from sympy.combinatorics.perm_groups import PermutationGroup
-        >>> from sympy.combinatorics.util import _verify_bsgs
+        >>> from sympy.combinatorics.testutil import _verify_bsgs
         >>> S = SymmetricGroup(7)
         >>> prop_even = lambda x: x.is_even
         >>> base, strong_gens = S.schreier_sims_incremental()
@@ -2585,7 +3054,7 @@ def subgroup_search(self, prop, base=None, strong_gens=None, tests=None,\
         # compute BSGS-related structures
         strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
         basic_orbits, transversals = _orbits_transversals_from_bsgs(base,\
-                                                                    strong_gens_distr)
+                                     strong_gens_distr)
         # handle subgroup initialization and tests
         if init_subgroup is None:
             init_subgroup = PermutationGroup([identity])
@@ -2602,7 +3071,8 @@ def subgroup_search(self, prop, base=None, strong_gens=None, tests=None,\
         res_base = base[:]
         # line 3: compute BSGS and related structures for K
         res_base, res_strong_gens = res.schreier_sims_incremental(base=res_base)
-        res_strong_gens_distr = _distribute_gens_by_base(res_base, res_strong_gens)
+        res_strong_gens_distr = _distribute_gens_by_base(res_base,\
+                                res_strong_gens)
         res_basic_orbits_init_base =\
         [PermutationGroup(res_strong_gens_distr[i]).orbit(res_base[i])\
          for i in range(base_len)]
@@ -2649,7 +3119,7 @@ def subgroup_search(self, prop, base=None, strong_gens=None, tests=None,\
                   base_ordering[mu[l]] and\
                   base_ordering[computed_words[l](base[l])] <\
                   base_ordering[nu[l]] and\
-                  tests[l](computed_words[base_len - 1]):
+                  tests[l](computed_words):
                 # line 11: change the (partial) base of K
                 new_point = computed_words[l](base[l])
                 res_base[l] = new_point
@@ -2699,7 +3169,7 @@ def subgroup_search(self, prop, base=None, strong_gens=None, tests=None,\
                base_ordering[temp_point] > base_ordering[mu[l]] and\
                base_ordering[temp_point] < base_ordering[nu[l]] and\
                temp_point in orbit_reps[l] and\
-               tests[l](g) and\
+               tests[l](computed_words) and\
                prop(g):
                 # line 19: reset the base of K
                 gens = res.generators[:]
@@ -2754,7 +3224,10 @@ def subgroup_search(self, prop, base=None, strong_gens=None, tests=None,\
             # line 29: set the next element from the current branch and update
             # accorndingly
             c[l] += 1
-            element = ~(computed_words[l - 1])
+            if l == 0:
+                element = identity
+            else:
+                element = ~(computed_words[l - 1])
             gamma  = element(sorted_orbits[l][c[l]])
             u[l] = transversals[l][gamma]
             if l == 0:
@@ -2805,55 +3278,4 @@ def transitivity_degree(self):
             return self._transitivity_degree
 
 
-def DirectProduct(*groups):
-    """
-    Returns the direct product of several groups as a permutation group.
-
-    This is implemented much like the __mul__ procedure for taking the direct
-    product of two permutation groups, but the idea of shifting the
-    generators is realized in the case of an arbitrary number of groups.
-    A call to DirectProduct(G1, G2, ..., Gn) is generally expected to be faster
-    than a call to G1*G2*...*Gn (and thus the need for this algorithm).
-
-    Examples
-    ========
-
-    >>> from sympy.combinatorics.perm_groups import DirectProduct
-    >>> from sympy.combinatorics.named_groups import CyclicGroup
-    >>> C = CyclicGroup(4)
-    >>> G = DirectProduct(C,C,C)
-    >>> G.order()
-    64
-
-    See Also
-    ========
-    __mul__
-
-    """
-    degrees = []
-    gens_count = []
-    total_degree = 0
-    total_gens = 0
-    for group in groups:
-        current_deg = group.degree
-        current_num_gens = len(group.generators)
-        degrees.append(current_deg)
-        total_degree += current_deg
-        gens_count.append(current_num_gens)
-        total_gens += current_num_gens
-    array_gens = []
-    for i in range(total_gens):
-        array_gens.append(range(total_degree))
-    current_gen = 0
-    current_deg = 0
-    for i in xrange(len(gens_count)):
-        for j in xrange(current_gen, current_gen + gens_count[i]):
-            gen = ((groups[i].generators)[j - current_gen]).array_form
-            array_gens[j][current_deg:current_deg + degrees[i]] =\
-            [ x + current_deg for x in gen]
-        current_gen += gens_count[i]
-        current_deg += degrees[i]
-    perm_gens = [_new_from_array_form(array) for array in array_gens]
-    return PermutationGroup(perm_gens)
-
 PermGroup = PermutationGroup
diff --git a/sympy/combinatorics/tests/test_group_constructs.py b/sympy/combinatorics/tests/test_group_constructs.py
new file mode 100644
index 000000000000..f033c47ab808
--- /dev/null
+++ b/sympy/combinatorics/tests/test_group_constructs.py
@@ -0,0 +1,14 @@
+from sympy.combinatorics.group_constructs import DirectProduct
+from sympy.combinatorics.named_groups import CyclicGroup, DihedralGroup
+
+def test_direct_product_n():
+    C = CyclicGroup(4)
+    D = DihedralGroup(4)
+    G = DirectProduct(C, C, C)
+    assert G.order() == 64
+    assert G.degree == 12
+    assert len(G.orbits()) == 3
+    assert G.is_abelian == True
+    H = DirectProduct(D, C)
+    assert H.order() == 32
+    assert H.is_abelian == False
diff --git a/sympy/combinatorics/tests/test_named_groups.py b/sympy/combinatorics/tests/test_named_groups.py
index 9adb32569fe7..fb2ce28e6cae 100644
--- a/sympy/combinatorics/tests/test_named_groups.py
+++ b/sympy/combinatorics/tests/test_named_groups.py
@@ -7,7 +7,7 @@ def test_SymmetricGroup():
     elements = list(G.generate())
     assert (G.generators[0]).size == 5
     assert len(elements) == 120
-    assert G.is_solvable() == False
+    assert G.is_solvable == False
     assert G.is_abelian == False
     assert G.is_transitive == True
     H = SymmetricGroup(1)
@@ -19,7 +19,7 @@ def test_CyclicGroup():
     G = CyclicGroup(10)
     elements = list(G.generate())
     assert len(elements) == 10
-    assert (G.commutator()).order() == 1
+    assert (G.derived_subgroup()).order() == 1
     assert G.is_abelian == True
     H = CyclicGroup(1)
     assert H.order() == 1
diff --git a/sympy/combinatorics/tests/test_perm_groups.py b/sympy/combinatorics/tests/test_perm_groups.py
index 52883abadaec..0bd77293a454 100644
--- a/sympy/combinatorics/tests/test_perm_groups.py
+++ b/sympy/combinatorics/tests/test_perm_groups.py
@@ -1,12 +1,14 @@
-from sympy.combinatorics.perm_groups import PermutationGroup, DirectProduct
+from sympy.combinatorics.perm_groups import PermutationGroup
+from sympy.combinatorics.group_constructs import DirectProduct
 from sympy.combinatorics.named_groups import SymmetricGroup, CyclicGroup,\
-DihedralGroup, AlternatingGroup
+DihedralGroup, AlternatingGroup, AbelianGroup
 from sympy.combinatorics.permutations import Permutation, perm_af_muln, cyclic
 from sympy.utilities.pytest import raises, skip, XFAIL
 from sympy.combinatorics.generators import rubik_cube_generators
 import random
-from sympy.combinatorics.util import _verify_bsgs
-
+from sympy.combinatorics.testutil import _verify_bsgs, _verify_centralizer,\
+_cmp_perm_lists, _verify_normal_closure
+from sympy.combinatorics.util import _distribute_gens_by_base
 
 def test_new():
     a = Permutation([1, 0])
@@ -29,7 +31,6 @@ def test_generate():
     assert list(g) == [Permutation([0, 1]), Permutation([1, 0])]
     g = PermutationGroup([a]).generate(method='dimino')
     assert list(g) == [Permutation([0, 1]), Permutation([1, 0])]
-
     a = Permutation([2, 0, 1])
     b = Permutation([2, 1, 0])
     G = PermutationGroup([a, b])
@@ -37,10 +38,8 @@ def test_generate():
     v1 = [p.array_form for p in list(g)]
     v1.sort()
     assert v1 == [[0,1,2], [0,2,1], [1,0,2], [1,2,0], [2,0,1], [2,1,0]]
-
     v2 = list(G.generate(method='dimino', af=True))
     assert v1 == sorted(v2)
-
     a = Permutation([2, 0, 1, 3, 4, 5])
     b = Permutation([2, 1, 3, 4, 5, 0])
     g = PermutationGroup([a, b]).generate(af=True)
@@ -76,6 +75,77 @@ def test_stabilizer():
     G2 = G.stabilizer(2)
     assert G2.order() == 181440
 
+def test_center():
+    # the center of the dihedral group D_n is of order 2 for even n
+    for i in (4, 6, 10):
+        D = DihedralGroup(i)
+        assert (D.center()).order() == 2
+    # the center of the dihedral group D_n is of order 1 for odd n>2
+    for i in (3, 5, 7):
+        D = DihedralGroup(i)
+        assert (D.center()).order() == 1
+    # the center of an abelian group is the group itself
+    for i in (2, 3, 5):
+        for j in (1, 5, 7):
+            for k in (1, 1, 11):
+                G = AbelianGroup(i, j, k)
+                assert G.center() == G
+    # the center of a nonabelian simple group is trivial
+    for i in(1, 5, 9):
+        A = AlternatingGroup(i)
+        assert (A.center()).order() == 1
+    # brute-force verifications
+    D = DihedralGroup(5)
+    A = AlternatingGroup(3)
+    C = CyclicGroup(4)
+    G = D*A*C
+    assert _verify_centralizer(G, G)
+
+def test_centralizer():
+    # the centralizer of the trivial group is the entire group
+    S = SymmetricGroup(2)
+    assert S.centralizer(Permutation(range(2))) == S
+    A = AlternatingGroup(5)
+    assert A.centralizer(Permutation(range(5))) == A
+    # a centralizer in the trivial group is the trivial group itself
+    triv = PermutationGroup([Permutation([0,1,2,3])])
+    D = DihedralGroup(4)
+    assert triv.centralizer(D) == triv
+    # brute-force verifications for centralizers of groups
+    for i in (4, 5, 6):
+        S = SymmetricGroup(i)
+        A = AlternatingGroup(i)
+        C = CyclicGroup(i)
+        D = DihedralGroup(i)
+        for gp in (S, A, C, D):
+            for gp2 in (S, A, C, D):
+                if gp2 != gp:
+                    assert _verify_centralizer(gp, gp2)
+    # verify the centralizer for all elements of several groups
+    S = SymmetricGroup(5)
+    elements = list(S.generate_dimino())
+    for element in elements:
+        assert _verify_centralizer(S, element)
+    A = AlternatingGroup(5)
+    elements = list(A.generate_dimino())
+    for element in elements:
+        assert _verify_centralizer(A, element)
+    D = DihedralGroup(7)
+    elements = list(D.generate_dimino())
+    for element in elements:
+        assert _verify_centralizer(D, element)
+    # verify centralizers of small groups within small groups
+    small = []
+    for i in (1, 2, 3):
+        small.append(SymmetricGroup(i))
+        small.append(AlternatingGroup(i))
+        small.append(DihedralGroup(i))
+        small.append(CyclicGroup(i))
+    for gp in small:
+        for gp2 in small:
+            if gp.degree == gp2.degree:
+                assert _verify_centralizer(gp, gp2)
+
 def test_coset_repr():
     a = Permutation([0, 2, 1])
     b = Permutation([1, 0, 2])
@@ -175,11 +245,11 @@ def test_eq():
     assert G1 != G4
     assert not G4.is_subgroup(G1)
 
-def test_commutators():
+def test_derived_subgroup():
     a = Permutation([1, 0, 2, 4, 3])
     b = Permutation([0, 1, 3, 2, 4])
     G = PermutationGroup([a,b])
-    C = G.commutator()
+    C = G.derived_subgroup()
     assert C.order() == 3
     assert C.is_normal(G)
     assert C.is_subgroup(G)
@@ -187,18 +257,18 @@ def test_commutators():
     gens_cube = [[1, 3, 5, 7, 0, 2, 4, 6], [1, 3, 0, 2, 5, 7, 4, 6]]
     gens = [Permutation(p) for p in gens_cube]
     G = PermutationGroup(gens)
-    C = G.commutator()
+    C = G.derived_subgroup()
     assert C.order() == 12
 
 def test_is_solvable():
     a = Permutation([1,2,0])
     b = Permutation([1,0,2])
     G = PermutationGroup([a, b])
-    assert G.is_solvable()
+    assert G.is_solvable
     a = Permutation([1,2,3,4,0])
     b = Permutation([1,0,2,3,4])
     G = PermutationGroup([a, b])
-    assert not G.is_solvable()
+    assert not G.is_solvable
 
 def test_rubik1():
     gens = rubik_cube_generators()
@@ -209,7 +279,7 @@ def test_rubik1():
     G2 = PermutationGroup(gens2)
     assert G2.order() == 663552
     assert G2.is_subgroup(G1)
-    C1 = G1.commutator()
+    C1 = G1.derived_subgroup()
     assert C1.order() == 4877107200
     assert C1.is_subgroup(G1)
     assert not G2.is_subgroup(C1)
@@ -365,18 +435,6 @@ def test_baseswap():
     assert randomized[0] == [0, 2, 1]
     assert _verify_bsgs(S, randomized[0], randomized[1]) == True
 
-def test_direct_product_n():
-    C = CyclicGroup(4)
-    D = DihedralGroup(4)
-    G = DirectProduct(C, C, C)
-    assert G.order() == 64
-    assert G.degree == 12
-    assert len(G.orbits()) == 3
-    assert G.is_abelian == True
-    H = DirectProduct(D, C)
-    assert H.order() == 32
-    assert H.is_abelian == False
-
 def test_schreier_sims_incremental():
     identity = Permutation([0, 1, 2, 3, 4])
     TrivialGroup = PermutationGroup([identity])
@@ -420,7 +478,8 @@ def test_subgroup_search():
         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)
+        assert S.stabilizer(3).stabilizer(4) ==\
+               S.subgroup_search(prop_fix_points)
         points = [3, 5]
         fix35 = A.subgroup_search(prop_fix_points)
         points = [5]
@@ -428,6 +487,123 @@ def test_subgroup_search():
         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)
+        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]
+
+def test_normal_closure():
+    # the normal closure of the trivial group is trivial
+    S = SymmetricGroup(3)
+    identity = Permutation([0, 1, 2])
+    closure = S.normal_closure(identity)
+    assert closure.is_trivial
+    # the normal closure of the entire group is the entire group
+    A = AlternatingGroup(4)
+    assert A.normal_closure(A) == A
+    # brute-force verifications for subgroups
+    for i in (3, 4, 5):
+        S = SymmetricGroup(i)
+        A = AlternatingGroup(i)
+        D = DihedralGroup(i)
+        C = CyclicGroup(i)
+        for gp in (A, D, C):
+            assert _verify_normal_closure(S, gp)
+    # brute-force verifications for all elements of a group
+    S = SymmetricGroup(5)
+    elements = list(S.generate_dimino())
+    for element in elements:
+        assert _verify_normal_closure(S, element)
+    # small groups
+    small = []
+    for i in (1, 2, 3):
+        small.append(SymmetricGroup(i))
+        small.append(AlternatingGroup(i))
+        small.append(DihedralGroup(i))
+        small.append(CyclicGroup(i))
+    for gp in small:
+        for gp2 in small:
+            if gp2.is_subgroup(gp):
+                assert _verify_normal_closure(gp, gp2)
+
+def test_derived_series():
+    # the derived series of the trivial group consists only of the trivial group
+    triv = PermutationGroup([Permutation([0, 1, 2])])
+    assert triv.derived_series() == [triv]
+    # the derived series for a simple group consists only of the group itself
+    for i in (5, 6, 7):
+        A = AlternatingGroup(i)
+        assert A.derived_series() == [A]
+    # the derived series for S_4 is S_4 > A_4 > K_4 > triv
+    S = SymmetricGroup(4)
+    series = S.derived_series()
+    assert series[1] == AlternatingGroup(4)
+    assert series[2] == DihedralGroup(2)
+    assert series[3].is_trivial
+
+def test_lower_central_series():
+    # the lower central series of the trivial group consists of the trivial
+    # group
+    triv = PermutationGroup([Permutation([0, 1, 2])])
+    assert triv.lower_central_series() == [triv]
+    # the lower central series of a simple group consists of the group itself
+    for i in (5, 6, 7):
+        A = AlternatingGroup(i)
+        assert A.lower_central_series() == [A]
+    # GAP-verified example
+    S = SymmetricGroup(6)
+    series = S.lower_central_series()
+    assert len(series) == 2
+    assert series[1] == AlternatingGroup(6)
+
+def test_commutator():
+    # the commutator of the trivial group and the trivial group is trivial
+    S = SymmetricGroup(3)
+    triv = PermutationGroup([Permutation([0, 1, 2])])
+    assert S.commutator(triv, triv) == triv
+    # the commutator of the trivial group and any other group is again trivial
+    A = AlternatingGroup(3)
+    assert S.commutator(triv, A) == triv
+    # the commutator is commutative
+    for i in (3, 4, 5):
+        S = SymmetricGroup(i)
+        A = AlternatingGroup(i)
+        D = DihedralGroup(i)
+        assert S.commutator(A, D) == S.commutator(D, A)
+    # the commutator of an abelian group is trivial
+    S = SymmetricGroup(7)
+    A1 = AbelianGroup(2, 5)
+    A2 = AbelianGroup(3, 4)
+    triv = PermutationGroup([Permutation([0, 1, 2, 3, 4, 5, 6])])
+    assert S.commutator(A1, A1) == triv
+    assert S.commutator(A2, A2) == triv
+    # examples calculated by hand
+    S = SymmetricGroup(3)
+    A = AlternatingGroup(3)
+    assert S.commutator(A, S) == A
+
+def test_is_nilpotent():
+    # every abelian group is nilpotent
+    for i in (1, 2, 3):
+        C = CyclicGroup(i)
+        Ab = AbelianGroup(i, i + 2)
+        assert C.is_nilpotent
+        assert Ab.is_nilpotent
+    Ab = AbelianGroup(5, 7, 10)
+    assert Ab.is_nilpotent
+    # A_5 is not solvable and thus not nilpotent
+    assert AlternatingGroup(5).is_nilpotent == False
+
+def test_is_trivial():
+    for i in range(5):
+        triv = PermutationGroup([Permutation(range(i))])
+        assert triv.is_trivial
+
+def test_pointwise_stabilizer():
+    S = SymmetricGroup(5)
+    points = []
+    stab = S
+    for point in (2, 0, 3, 4, 1):
+        stab = stab.stabilizer(point)
+        points.append(point)
+        assert S.pointwise_stabilizer(points) == stab
diff --git a/sympy/combinatorics/tests/test_testutil.py b/sympy/combinatorics/tests/test_testutil.py
new file mode 100644
index 000000000000..63a1a1bdb1ac
--- /dev/null
+++ b/sympy/combinatorics/tests/test_testutil.py
@@ -0,0 +1,52 @@
+from sympy.combinatorics.named_groups import SymmetricGroup, AlternatingGroup,\
+CyclicGroup
+from sympy.combinatorics.testutil import _verify_bsgs, _cmp_perm_lists,\
+_naive_list_centralizer, _verify_centralizer,\
+_verify_normal_closure
+from sympy.combinatorics.permutations import Permutation
+from sympy.combinatorics.perm_groups import PermutationGroup
+from random import shuffle
+
+
+def test_cmp_perm_lists():
+    S = SymmetricGroup(4)
+    els = list(S.generate_dimino())
+    other = els[:]
+    shuffle(other)
+    assert _cmp_perm_lists(els, other) == True
+
+def test_naive_list_centralizer():
+    # verified by GAP
+    S = SymmetricGroup(3)
+    A = AlternatingGroup(3)
+    assert _naive_list_centralizer(S, S) == [Permutation([0, 1, 2])]
+    assert PermutationGroup(_naive_list_centralizer(S, A)) == A
+
+def test_verify_bsgs():
+    S = SymmetricGroup(5)
+    S.schreier_sims()
+    base = S.base
+    strong_gens = S.strong_gens
+    gens = S.generators
+    assert _verify_bsgs(S, base, strong_gens) == True
+    assert _verify_bsgs(S, base[:-1], strong_gens) == False
+    assert _verify_bsgs(S, base, S.generators) == False
+
+def test_verify_centralizer():
+    # verified by GAP
+    S = SymmetricGroup(3)
+    A = AlternatingGroup(3)
+    triv = PermutationGroup([Permutation([0, 1, 2])])
+    assert _verify_centralizer(S, S, centr=triv)
+    assert _verify_centralizer(S, A, centr=A)
+
+def test_verify_normal_closure():
+    # verified by GAP
+    S = SymmetricGroup(3)
+    A = AlternatingGroup(3)
+    assert _verify_normal_closure(S, A, closure=A)
+    S = SymmetricGroup(5)
+    A = AlternatingGroup(5)
+    C = CyclicGroup(5)
+    assert _verify_normal_closure(S, A, closure=A)
+    assert _verify_normal_closure(S, C, closure=A)
diff --git a/sympy/combinatorics/tests/test_util.py b/sympy/combinatorics/tests/test_util.py
index 72558e8925a8..f6c780ddf70d 100644
--- a/sympy/combinatorics/tests/test_util.py
+++ b/sympy/combinatorics/tests/test_util.py
@@ -5,7 +5,8 @@
 from sympy.combinatorics.util import _check_cycles_alt_sym, _strip,\
 _distribute_gens_by_base, _strong_gens_from_distr,\
 _orbits_transversals_from_bsgs, _handle_precomputed_bsgs, _base_ordering,\
-_verify_bsgs, _remove_gens
+_remove_gens
+from sympy.combinatorics.testutil import _verify_bsgs
 
 def test_check_cycles_alt_sym():
     perm1 = Permutation([[0, 1, 2, 3, 4, 5, 6], [7], [8], [9]])
@@ -97,16 +98,6 @@ def test_base_ordering():
     degree = 7
     assert _base_ordering(base, degree) == [3, 4, 0, 5, 1, 2, 6]
 
-def test_verify_bsgs():
-    S = SymmetricGroup(5)
-    S.schreier_sims()
-    base = S.base
-    strong_gens = S.strong_gens
-    gens = S.generators
-    assert _verify_bsgs(S, base, strong_gens) == True
-    assert _verify_bsgs(S, base[:-1], strong_gens) == False
-    assert _verify_bsgs(S, base, S.generators) == False
-
 def test_remove_gens():
     S = SymmetricGroup(10)
     base, strong_gens = S.schreier_sims_incremental()
diff --git a/sympy/combinatorics/testutil.py b/sympy/combinatorics/testutil.py
new file mode 100644
index 000000000000..afa4a8848b2f
--- /dev/null
+++ b/sympy/combinatorics/testutil.py
@@ -0,0 +1,184 @@
+from sympy.combinatorics.named_groups import SymmetricGroup, DihedralGroup,\
+AlternatingGroup, CyclicGroup
+from sympy.combinatorics.util import _distribute_gens_by_base
+
+def _cmp_perm_lists(first, second):
+    """
+    Compare two lists of permutations as sets.
+
+    This is used for testing purposes. Since the array form of a
+    permutation is currently a list, Permutation is not hashable
+    and cannot be put into a set.
+
+    Examples
+    ========
+
+    >>> from sympy.combinatorics.permutations import Permutation
+    >>> from sympy.combinatorics.testutil import _cmp_perm_lists
+    >>> a = Permutation([0, 2, 3, 4, 1])
+    >>> b = Permutation([1, 2, 0, 4, 3])
+    >>> c = Permutation([3, 4, 0, 1, 2])
+    >>> ls1 = [a, b, c]
+    >>> ls2 = [b, c, a]
+    >>> _cmp_perm_lists(ls1, ls2)
+    True
+
+    """
+    first.sort(key = lambda x: x.array_form)
+    second.sort(key = lambda x: x.array_form)
+    return first == second
+
+def _naive_list_centralizer(self, other):
+    from sympy.combinatorics.perm_groups import PermutationGroup
+    """
+    Return a list of elements for the centralizer of a subgroup/set/element.
+
+    This is a brute-force implementation that goes over all elements of the
+    group and checks for membership in the centralizer. It is used to
+    test ``.centralizer()`` from ``sympy.combinatorics.perm_groups``.
+
+    Examples
+    ========
+    >>> from sympy.combinatorics.testutil import _naive_list_centralizer
+    >>> from sympy.combinatorics.named_groups import DihedralGroup
+    >>> D = DihedralGroup(4)
+    >>> _naive_list_centralizer(D, D)
+    [Permutation([0, 1, 2, 3]), Permutation([2, 3, 0, 1])]
+
+    See Also
+    ========
+
+    sympy.combinatorics.perm_groups.centralizer
+
+    """
+    if hasattr(other, 'generators'):
+        elements = list(self.generate_dimino())
+        gens = other.generators
+        commutes_with_gens = lambda x: [x*gen for gen in gens] ==\
+                                       [gen*x for gen in gens]
+        centralizer_list = []
+        for element in elements:
+            if commutes_with_gens(element):
+                centralizer_list.append(element)
+        return centralizer_list
+    elif hasattr(other, 'getitem'):
+        return _naive_list_centralizer(self, PermutationGroup(other))
+    elif hasattr(other, 'array_form'):
+        return _naive_list_centralizer(self, PermutationGroup([other]))
+
+def _verify_bsgs(group, base, gens):
+    """
+    Verify the correctness of a base and strong generating set.
+
+    This is a naive implementation using the definition of a base and a strong
+    generating set relative to it. There are other procedures for
+    verifying a base and strong generating set, but this one will
+    serve for more robust testing.
+
+    Examples
+    ========
+
+    >>> from sympy.combinatorics.named_groups import AlternatingGroup
+    >>> from sympy.combinatorics.testutil import _verify_bsgs
+    >>> A = AlternatingGroup(4)
+    >>> A.schreier_sims()
+    >>> _verify_bsgs(A, A.base, A.strong_gens)
+    True
+
+    See Also
+    ========
+
+    sympy.combinatorics.perm_groups.PermutationGroup.schreier_sims
+
+    """
+    from sympy.combinatorics.perm_groups import PermutationGroup
+    strong_gens_distr = _distribute_gens_by_base(base, gens)
+    base_len = len(base)
+    degree = group.degree
+    current_stabilizer = group
+    for i in range(base_len):
+        candidate = PermutationGroup(strong_gens_distr[i])
+        if current_stabilizer.order() != candidate.order():
+            return False
+        current_stabilizer = current_stabilizer.stabilizer(base[i])
+    if current_stabilizer.order() != 1:
+        return False
+    return True
+
+def _verify_centralizer(group, arg, centr=None):
+    """
+    Verify the centralizer of a group/set/element inside another group.
+
+    This is used for testing ``.centralizer()`` from
+    ``sympy.combinatorics.perm_groups``
+
+    Examples
+    ========
+
+    >>> from sympy.combinatorics.named_groups import (SymmetricGroup,
+    ... AlternatingGroup)
+    >>> from sympy.combinatorics.perm_groups import PermutationGroup
+    >>> from sympy.combinatorics.permutations import Permutation
+    >>> from sympy.combinatorics.testutil import _verify_centralizer
+    >>> S = SymmetricGroup(5)
+    >>> A = AlternatingGroup(5)
+    >>> centr = PermutationGroup([Permutation([0, 1, 2, 3, 4])])
+    >>> _verify_centralizer(S, A, centr)
+    True
+
+    See Also
+    ========
+
+    _naive_list_centralizer,\
+    sympy.combinatorics.perm_groups.PermutationGroup.centralizer,\
+    _cmp_perm_lists
+
+    """
+    if centr is None:
+        centr = group.centralizer(arg)
+    centr_list = list(centr.generate_dimino())
+    centr_list_naive = _naive_list_centralizer(group, arg)
+    return _cmp_perm_lists(centr_list, centr_list_naive)
+
+def _verify_normal_closure(group, arg, closure=None):
+    from sympy.combinatorics.perm_groups import PermutationGroup
+    """
+    Verify the normal closure of a subgroup/subset/element in a group.
+
+    This is used to test
+    sympy.combinatorics.perm_groups.PermutationGroup.normal_closure
+
+    Examples
+    ========
+
+    >>> from sympy.combinatorics.named_groups import (SymmetricGroup,\
+    ... AlternatingGroup)
+    >>> from sympy.combinatorics.testutil import _verify_normal_closure
+    >>> S = SymmetricGroup(3)
+    >>> A = AlternatingGroup(3)
+    >>> _verify_normal_closure(S, A, closure=A)
+    True
+
+    See Also
+    ========
+
+    sympy.combinatorics.perm_groups.PermutationGroup.normal_closure
+
+    """
+    if closure is None:
+        closure = group.normal_closure(arg)
+    conjugates = []
+    group_els = list(group.generate_dimino())
+    if hasattr(arg, 'generators'):
+        subgr_gens = arg.generators
+    elif hasattr(arg, '__getitem__'):
+        subgr_gens = arg
+    elif hasattr(arg, 'array_form'):
+        subgr_gens = [arg]
+    for el in group_els:
+        for gen in subgr_gens:
+            conjugate = (~el)*gen*el
+            if conjugate not in conjugates:
+                conjugates.append(conjugate)
+    naive_closure = PermutationGroup(conjugates)
+    return closure == naive_closure
diff --git a/sympy/combinatorics/util.py b/sympy/combinatorics/util.py
index 4ac249631c90..4370495a8b7a 100644
--- a/sympy/combinatorics/util.py
+++ b/sympy/combinatorics/util.py
@@ -193,10 +193,10 @@ def _handle_precomputed_bsgs(base, strong_gens, transversals=None,\
     Returns
     =======
 
-    ``(transversals, basic_orbits, strong_gens_distr)`` where ``transversals`` are the
-    basic transversals, ``basic_orbits`` are the basic orbits, and
-    ``strong_gens_distr`` are the strong generators distributed by membership in basic
-    stabilizers.
+    ``(transversals, basic_orbits, strong_gens_distr)`` where ``transversals``
+    are the basic transversals, ``basic_orbits`` are the basic orbits, and
+    ``strong_gens_distr`` are the strong generators distributed by membership
+    in basic stabilizers.
 
     Examples
     ========
@@ -314,7 +314,8 @@ def _remove_gens(base, strong_gens, basic_orbits=None, strong_gens_distr=None):
 
     >>> from sympy.combinatorics.named_groups import SymmetricGroup
     >>> from sympy.combinatorics.perm_groups import PermutationGroup
-    >>> from sympy.combinatorics.util import _remove_gens, _verify_bsgs
+    >>> from sympy.combinatorics.util import _remove_gens
+    >>> from sympy.combinatorics.testutil import _verify_bsgs
     >>> S = SymmetricGroup(15)
     >>> base, strong_gens = S.schreier_sims_incremental()
     >>> len(strong_gens)
@@ -481,42 +482,3 @@ def _strong_gens_from_distr(strong_gens_distr):
             if gen not in result:
                 result.append(gen)
         return result
-
-def _verify_bsgs(group, base, gens):
-    """
-    Verify the correctness of a base and strong generating set.
-
-    This is a naive implementation using the definition of a base and a strong
-    generating set relative to it. There are other procedures for
-    verifying a base and strong generating set, but this one will
-    serve for more robust testing.
-
-    Examples
-    ========
-
-    >>> from sympy.combinatorics.named_groups import AlternatingGroup
-    >>> from sympy.combinatorics.util import _verify_bsgs
-    >>> A = AlternatingGroup(4)
-    >>> A.schreier_sims()
-    >>> _verify_bsgs(A, A.base, A.strong_gens)
-    True
-
-    See Also
-    ========
-
-    sympy.combinatorics.perm_groups.PermutationGroup.schreier_sims
-
-    """
-    from sympy.combinatorics.perm_groups import PermutationGroup
-    strong_gens_distr = _distribute_gens_by_base(base, gens)
-    base_len = len(base)
-    degree = group.degree
-    current_stabilizer = group
-    for i in range(base_len):
-        candidate = PermutationGroup(strong_gens_distr[i])
-        if current_stabilizer.order() != candidate.order():
-            return False
-        current_stabilizer = current_stabilizer.stabilizer(base[i])
-    if current_stabilizer.order() != 1:
-        return False
-    return True