text: from __future__ import print_function, division
  
  from random import randrange, choice
  from math import log
  
  from sympy.core import Basic
  from sympy.core.compatibility import range
  from sympy.combinatorics import Permutation
  from sympy.combinatorics.permutations import (_af_commutes_with, _af_invert,
      _af_rmul, _af_rmuln, _af_pow, Cycle)
  from sympy.combinatorics.util import (_check_cycles_alt_sym,
      _distribute_gens_by_base, _orbits_transversals_from_bsgs,
      _handle_precomputed_bsgs, _base_ordering, _strong_gens_from_distr,
      _strip, _strip_af)
  from sympy.functions.combinatorial.factorials import factorial
  from sympy.ntheory import sieve
  from sympy.utilities.iterables import has_variety, is_sequence, uniq
  from sympy.utilities.randtest import _randrange
  from itertools import islice
  
  rmul = Permutation.rmul_with_af
  _af_new = Permutation._af_new
  
  
  class PermutationGroup(Basic):
      """The class defining a Permutation group.
  
      PermutationGroup([p1, p2, ..., pn]) returns the permutation group
      generated by the list of permutations. This group can be supplied
      to Polyhedron if one desires to decorate the elements to which the
      indices of the permutation refer.
  
      Examples
      ========
  
      >>> from sympy.combinatorics import Permutation
      >>> Permutation.print_cyclic = True
      >>> from sympy.combinatorics.permutations import Cycle
      >>> from sympy.combinatorics.polyhedron import Polyhedron
      >>> from sympy.combinatorics.perm_groups import PermutationGroup
  
      The permutations corresponding to motion of the front, right and
      bottom face of a 2x2 Rubik's cube are defined:
  
      >>> F = Permutation(2, 19, 21, 8)(3, 17, 20, 10)(4, 6, 7, 5)
      >>> R = Permutation(1, 5, 21, 14)(3, 7, 23, 12)(8, 10, 11, 9)
      >>> D = Permutation(6, 18, 14, 10)(7, 19, 15, 11)(20, 22, 23, 21)
  
      These are passed as permutations to PermutationGroup:
  
      >>> G = PermutationGroup(F, R, D)
      >>> G.order()
      3674160
  
      The group can be supplied to a Polyhedron in order to track the
      objects being moved. An example involving the 2x2 Rubik's cube is
      given there, but here is a simple demonstration:
  
      >>> a = Permutation(2, 1)
      >>> b = Permutation(1, 0)
      >>> G = PermutationGroup(a, b)
      >>> P = Polyhedron(list('ABC'), pgroup=G)
      >>> P.corners
      (A, B, C)
      >>> P.rotate(0) # apply permutation 0
      >>> P.corners
      (A, C, B)
      >>> P.reset()
      >>> P.corners
      (A, B, C)
  
      Or one can make a permutation as a product of selected permutations
      and apply them to an iterable directly:
  
      >>> P10 = G.make_perm([0, 1])
      >>> P10('ABC')
      ['C', 'A', 'B']
  
      See Also
      ========
  
      sympy.combinatorics.polyhedron.Polyhedron,
      sympy.combinatorics.permutations.Permutation
  
      References
      ==========
  
      [1] Holt, D., Eick, B., O'Brien, E.
      "Handbook of Computational Group Theory"
  
      [2] Seress, A.
      "Permutation Group Algorithms"
  
      [3] http://en.wikipedia.org/wiki/Schreier_vector
  
      [4] http://en.wikipedia.org/wiki/Nielsen_transformation
      #Product_replacement_algorithm
  
      [5] Frank Celler, Charles R.Leedham-Green, Scott H.Murray,
      Alice C.Niemeyer, and E.A.O'Brien. "Generating Random
      Elements of a Finite Group"
  
      [6] http://en.wikipedia.org/wiki/Block_%28permutation_group_theory%29
  
      [7] http://www.algorithmist.com/index.php/Union_Find
  
      [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
  
      [13] http://www.math.colostate.edu/~hulpke/CGT/cgtnotes.pdf
  
      """
      is_group = True
  
      def __new__(cls, *args, **kwargs):
          """The default constructor. Accepts Cycle and Permutation forms.
          Removes duplicates unless ``dups`` keyword is ``False``.
          """
          if not args:
              args = [Permutation()]
          else:
              args = list(args[0] if is_sequence(args[0]) else args)
          if any(isinstance(a, Cycle) for a in args):
              args = [Permutation(a) for a in args]
          if has_variety(a.size for a in args):
              degree = kwargs.pop('degree', None)
              if degree is None:
                  degree = max(a.size for a in args)
              for i in range(len(args)):
                  if args[i].size != degree:
                      args[i] = Permutation(args[i], size=degree)
          if kwargs.pop('dups', True):
              args = list(uniq([_af_new(list(a)) for a in args]))
          obj = Basic.__new__(cls, *args, **kwargs)
          obj._generators = args
          obj._order = None
          obj._center = []
          obj._is_abelian = None
          obj._is_transitive = None
          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
          obj._r = len(obj._generators)
          obj._degree = obj._generators[0].size
  
          # these attributes are assigned after running schreier_sims
          obj._base = []
          obj._strong_gens = []
          obj._strong_gens_slp = []
          obj._basic_orbits = []
          obj._transversals = []
          obj._transversal_slp = []
  
          # these attributes are assigned after running _random_pr_init
          obj._random_gens = []
          return obj
  
      def __getitem__(self, i):
          return self._generators[i]
  
      def __contains__(self, i):
          """Return ``True`` if `i` is contained in PermutationGroup.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation, PermutationGroup
          >>> p = Permutation(1, 2, 3)
          >>> Permutation(3) in PermutationGroup(p)
          True
  
          """
          if not isinstance(i, Permutation):
              raise TypeError("A PermutationGroup contains only Permutations as "
                              "elements, not elements of type %s" % type(i))
          return self.contains(i)
  
      def __len__(self):
          return len(self._generators)
  
      def __eq__(self, other):
          """Return ``True`` if PermutationGroup generated by elements in the
          group are same i.e they represent the same PermutationGroup.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> p = Permutation(0, 1, 2, 3, 4, 5)
          >>> G = PermutationGroup([p, p**2])
          >>> H = PermutationGroup([p**2, p])
          >>> G.generators == H.generators
          False
          >>> G == H
          True
  
          """
          if not isinstance(other, PermutationGroup):
              return False
  
          set_self_gens = set(self.generators)
          set_other_gens = set(other.generators)
  
          # before reaching the general case there are also certain
          # optimisation and obvious cases requiring less or no actual
          # computation.
          if set_self_gens == set_other_gens:
              return True
  
          # in the most general case it will check that each generator of
          # one group belongs to the other PermutationGroup and vice-versa
          for gen1 in set_self_gens:
              if not other.contains(gen1):
                  return False
          for gen2 in set_other_gens:
              if not self.contains(gen2):
                  return False
          return True
  
      def __hash__(self):
          return super(PermutationGroup, self).__hash__()
  
      def __mul__(self, other):
          """Return the direct product of two permutation groups as a permutation
          group.
  
          This implementation realizes the direct product by shifting the index
          set for the generators of the second group: so if we have `G` acting
          on `n1` points and `H` acting on `n2` points, `G*H` acts on `n1 + n2`
          points.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> from sympy.combinatorics.named_groups import CyclicGroup
          >>> G = CyclicGroup(5)
          >>> H = G*G
          >>> H
          PermutationGroup([
              (9)(0 1 2 3 4),
              (5 6 7 8 9)])
          >>> H.order()
          25
  
          """
          gens1 = [perm._array_form for perm in self.generators]
          gens2 = [perm._array_form for perm in other.generators]
          n1 = self._degree
          n2 = other._degree
          start = list(range(n1))
          end = list(range(n1, n1 + n2))
          for i in range(len(gens2)):
              gens2[i] = [x + n1 for x in gens2[i]]
          gens2 = [start + gen for gen in gens2]
          gens1 = [gen + end for gen in gens1]
          together = gens1 + gens2
          gens = [_af_new(x) for x in together]
          return PermutationGroup(gens)
  
      def _random_pr_init(self, r, n, _random_prec_n=None):
          r"""Initialize random generators for the product replacement algorithm.
  
          The implementation uses a modification of the original product
          replacement algorithm due to Leedham-Green, as described in [1],
          pp. 69-71; also, see [2], pp. 27-29 for a detailed theoretical
          analysis of the original product replacement algorithm, and [4].
  
          The product replacement algorithm is used for producing random,
          uniformly distributed elements of a group `G` with a set of generators
          `S`. For the initialization ``_random_pr_init``, a list ``R`` of
          `\max\{r, |S|\}` group generators is created as the attribute
          ``G._random_gens``, repeating elements of `S` if necessary, and the
          identity element of `G` is appended to ``R`` - we shall refer to this
          last element as the accumulator. Then the function ``random_pr()``
          is called ``n`` times, randomizing the list ``R`` while preserving
          the generation of `G` by ``R``. The function ``random_pr()`` itself
          takes two random elements ``g, h`` among all elements of ``R`` but
          the accumulator and replaces ``g`` with a randomly chosen element
          from `\{gh, g(~h), hg, (~h)g\}`. Then the accumulator is multiplied
          by whatever ``g`` was replaced by. The new value of the accumulator is
          then returned by ``random_pr()``.
  
          The elements returned will eventually (for ``n`` large enough) become
          uniformly distributed across `G` ([5]). For practical purposes however,
          the values ``n = 50, r = 11`` are suggested in [1].
  
          Notes
          =====
  
          THIS FUNCTION HAS SIDE EFFECTS: it changes the attribute
          self._random_gens
  
          See Also
          ========
  
          random_pr
  
          """
          deg = self.degree
          random_gens = [x._array_form for x in self.generators]
          k = len(random_gens)
          if k < r:
              for i in range(k, r):
                  random_gens.append(random_gens[i - k])
          acc = list(range(deg))
          random_gens.append(acc)
          self._random_gens = random_gens
  
          # handle randomized input for testing purposes
          if _random_prec_n is None:
              for i in range(n):
                  self.random_pr()
          else:
              for i in range(n):
                  self.random_pr(_random_prec=_random_prec_n[i])
  
      def _union_find_merge(self, first, second, ranks, parents, not_rep):
          """Merges two classes in a union-find data structure.
  
          Used in the implementation of Atkinson's algorithm as suggested in [1],
          pp. 83-87. The class merging process uses union by rank as an
          optimization. ([7])
  
          Notes
          =====
  
          THIS FUNCTION HAS SIDE EFFECTS: the list of class representatives,
          ``parents``, the list of class sizes, ``ranks``, and the list of
          elements that are not representatives, ``not_rep``, are changed due to
          class merging.
  
          See Also
          ========
  
          minimal_block, _union_find_rep
  
          References
          ==========
  
          [1] Holt, D., Eick, B., O'Brien, E.
          "Handbook of computational group theory"
  
          [7] http://www.algorithmist.com/index.php/Union_Find
  
          """
          rep_first = self._union_find_rep(first, parents)
          rep_second = self._union_find_rep(second, parents)
          if rep_first != rep_second:
              # union by rank
              if ranks[rep_first] >= ranks[rep_second]:
                  new_1, new_2 = rep_first, rep_second
              else:
                  new_1, new_2 = rep_second, rep_first
              total_rank = ranks[new_1] + ranks[new_2]
              if total_rank > self.max_div:
                  return -1
              parents[new_2] = new_1
              ranks[new_1] = total_rank
              not_rep.append(new_2)
              return 1
          return 0
  
      def _union_find_rep(self, num, parents):
          """Find representative of a class in a union-find data structure.
  
          Used in the implementation of Atkinson's algorithm as suggested in [1],
          pp. 83-87. After the representative of the class to which ``num``
          belongs is found, path compression is performed as an optimization
          ([7]).
  
          Notes
          =====
  
          THIS FUNCTION HAS SIDE EFFECTS: the list of class representatives,
          ``parents``, is altered due to path compression.
  
          See Also
          ========
  
          minimal_block, _union_find_merge
  
          References
          ==========
  
          [1] Holt, D., Eick, B., O'Brien, E.
          "Handbook of computational group theory"
  
          [7] http://www.algorithmist.com/index.php/Union_Find
  
          """
          rep, parent = num, parents[num]
          while parent != rep:
              rep = parent
              parent = parents[rep]
          # path compression
          temp, parent = num, parents[num]
          while parent != rep:
              parents[temp] = rep
              temp = parent
              parent = parents[temp]
          return rep
  
      @property
      def base(self):
          """Return a base from the Schreier-Sims algorithm.
  
          For a permutation group `G`, a base is a sequence of points
          `B = (b_1, b_2, ..., b_k)` such that no element of `G` apart
          from the identity fixes all the points in `B`. The concepts of
          a base and strong generating set and their applications are
          discussed in depth in [1], pp. 87-89 and [2], pp. 55-57.
  
          An alternative way to think of `B` is that it gives the
          indices of the stabilizer cosets that contain more than the
          identity permutation.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation, PermutationGroup
          >>> G = PermutationGroup([Permutation(0, 1, 3)(2, 4)])
          >>> G.base
          [0, 2]
  
          See Also
          ========
  
          strong_gens, basic_transversals, basic_orbits, basic_stabilizers
  
          """
          if self._base == []:
              self.schreier_sims()
          return self._base
  
      def baseswap(self, base, strong_gens, pos, randomized=False,
                   transversals=None, basic_orbits=None, strong_gens_distr=None):
          r"""Swap two consecutive base points in base and strong generating set.
  
          If a base for a group `G` is given by `(b_1, b_2, ..., b_k)`, this
          function returns a base `(b_1, b_2, ..., b_{i+1}, b_i, ..., b_k)`,
          where `i` is given by ``pos``, and a strong generating set relative
          to that base. The original base and strong generating set are not
          modified.
  
          The randomized version (default) is of Las Vegas type.
  
          Parameters
          ==========
  
          base, strong_gens
              The base and strong generating set.
          pos
              The position at which swapping is performed.
          randomized
              A switch between randomized and deterministic version.
          transversals
              The transversals for the basic orbits, if known.
          basic_orbits
              The basic orbits, if known.
          strong_gens_distr
              The strong generators distributed by basic stabilizers, if known.
  
          Returns
          =======
  
          (base, strong_gens)
              ``base`` is the new base, and ``strong_gens`` is a generating set
              relative to it.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.named_groups import SymmetricGroup
          >>> from sympy.combinatorics.testutil import _verify_bsgs
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> S = SymmetricGroup(4)
          >>> S.schreier_sims()
          >>> S.base
          [0, 1, 2]
          >>> base, gens = S.baseswap(S.base, S.strong_gens, 1, randomized=False)
          >>> base, gens
          ([0, 2, 1],
          [(0 1 2 3), (3)(0 1), (1 3 2),
           (2 3), (1 3)])
  
          check that base, gens is a BSGS
  
          >>> S1 = PermutationGroup(gens)
          >>> _verify_bsgs(S1, base, gens)
          True
  
          See Also
          ========
  
          schreier_sims
  
          Notes
          =====
  
          The deterministic version of the algorithm is discussed in
          [1], pp. 102-103; the randomized version is discussed in [1], p.103, and
          [2], p.98. It is of Las Vegas type.
          Notice that [1] contains a mistake in the pseudocode and
          discussion of BASESWAP: on line 3 of the pseudocode,
          `|\beta_{i+1}^{\left\langle T\right\rangle}|` should be replaced by
          `|\beta_{i}^{\left\langle T\right\rangle}|`, and the same for the
          discussion of the algorithm.
  
          """
          # construct the basic orbits, generators for the stabilizer chain
          # and transversal elements from whatever was provided
          transversals, basic_orbits, strong_gens_distr = \
              _handle_precomputed_bsgs(base, strong_gens, transversals,
                                   basic_orbits, strong_gens_distr)
          base_len = len(base)
          degree = self.degree
          # size of orbit of base[pos] under the stabilizer we seek to insert
          # in the stabilizer chain at position pos + 1
          size = len(basic_orbits[pos])*len(basic_orbits[pos + 1]) \
              //len(_orbit(degree, strong_gens_distr[pos], base[pos + 1]))
          # initialize the wanted stabilizer by a subgroup
          if pos + 2 > base_len - 1:
              T = []
          else:
              T = strong_gens_distr[pos + 2][:]
          # randomized version
          if randomized is True:
              stab_pos = PermutationGroup(strong_gens_distr[pos])
              schreier_vector = stab_pos.schreier_vector(base[pos + 1])
              # add random elements of the stabilizer until they generate it
              while len(_orbit(degree, T, base[pos])) != size:
                  new = stab_pos.random_stab(base[pos + 1],
                                             schreier_vector=schreier_vector)
                  T.append(new)
          # deterministic version
          else:
              Gamma = set(basic_orbits[pos])
              Gamma.remove(base[pos])
              if base[pos + 1] in Gamma:
                  Gamma.remove(base[pos + 1])
              # add elements of the stabilizer until they generate it by
              # ruling out member of the basic orbit of base[pos] along the way
              while len(_orbit(degree, T, base[pos])) != size:
                  gamma = next(iter(Gamma))
                  x = transversals[pos][gamma]
                  temp = x._array_form.index(base[pos + 1]) # (~x)(base[pos + 1])
                  if temp not in basic_orbits[pos + 1]:
                      Gamma = Gamma - _orbit(degree, T, gamma)
                  else:
                      y = transversals[pos + 1][temp]
                      el = rmul(x, y)
                      if el(base[pos]) not in _orbit(degree, T, base[pos]):
                          T.append(el)
                          Gamma = Gamma - _orbit(degree, T, base[pos])
          # build the new base and strong generating set
          strong_gens_new_distr = strong_gens_distr[:]
          strong_gens_new_distr[pos + 1] = T
          base_new = base[:]
          base_new[pos], base_new[pos + 1] = base_new[pos + 1], base_new[pos]
          strong_gens_new = _strong_gens_from_distr(strong_gens_new_distr)
          for gen in T:
              if gen not in strong_gens_new:
                  strong_gens_new.append(gen)
          return base_new, strong_gens_new
  
      @property
      def basic_orbits(self):
          """
          Return the basic orbits relative to a base and strong generating set.
  
          If `(b_1, b_2, ..., b_k)` is a base for a group `G`, and
          `G^{(i)} = G_{b_1, b_2, ..., b_{i-1}}` is the ``i``-th basic stabilizer
          (so that `G^{(1)} = G`), the ``i``-th basic orbit relative to this base
          is the orbit of `b_i` under `G^{(i)}`. See [1], pp. 87-89 for more
          information.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.named_groups import SymmetricGroup
          >>> S = SymmetricGroup(4)
          >>> S.basic_orbits
          [[0, 1, 2, 3], [1, 2, 3], [2, 3]]
  
          See Also
          ========
  
          base, strong_gens, basic_transversals, basic_stabilizers
  
          """
          if self._basic_orbits == []:
              self.schreier_sims()
          return self._basic_orbits
  
      @property
      def basic_stabilizers(self):
          """
          Return a chain of stabilizers relative to a base and strong generating
          set.
  
          The ``i``-th basic stabilizer `G^{(i)}` relative to a base
          `(b_1, b_2, ..., b_k)` is `G_{b_1, b_2, ..., b_{i-1}}`. For more
          information, see [1], pp. 87-89.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.named_groups import AlternatingGroup
          >>> A = AlternatingGroup(4)
          >>> A.schreier_sims()
          >>> A.base
          [0, 1]
          >>> for g in A.basic_stabilizers:
          ...     print(g)
          ...
          PermutationGroup([
              (3)(0 1 2),
              (1 2 3)])
          PermutationGroup([
              (1 2 3)])
  
          See Also
          ========
  
          base, strong_gens, basic_orbits, basic_transversals
  
          """
  
          if self._transversals == []:
              self.schreier_sims()
          strong_gens = self._strong_gens
          base = self._base
          if not base: # e.g. if self is trivial
              return []
          strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
          basic_stabilizers = []
          for gens in strong_gens_distr:
              basic_stabilizers.append(PermutationGroup(gens))
          return basic_stabilizers
  
      @property
      def basic_transversals(self):
          """
          Return basic transversals relative to a base and strong generating set.
  
          The basic transversals are transversals of the basic orbits. They
          are provided as a list of dictionaries, each dictionary having
          keys - the elements of one of the basic orbits, and values - the
          corresponding transversal elements. See [1], pp. 87-89 for more
          information.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.named_groups import AlternatingGroup
          >>> A = AlternatingGroup(4)
          >>> A.basic_transversals
          [{0: (3), 1: (3)(0 1 2), 2: (3)(0 2 1), 3: (0 3 1)}, {1: (3), 2: (1 2 3), 3: (1 3 2)}]
  
          See Also
          ========
  
          strong_gens, base, basic_orbits, basic_stabilizers
  
          """
  
          if self._transversals == []:
              self.schreier_sims()
          return self._transversals
  
      def coset_transversal(self, H):
          """Return a transversal of the right cosets of self by its subgroup H
          using the second method described in [1], Subsection 4.6.7
  
          """
  
          if not H.is_subgroup(self):
              raise ValueError("The argument must be a subgroup")
  
          if H.order() == 1:
              return self._elements
  
          self._schreier_sims(base=H.base) # make G.base an extension of H.base
  
          base = self.base
          base_ordering = _base_ordering(base, self.degree)
          identity = Permutation(self.degree - 1)
  
          transversals = self.basic_transversals[:]
          # transversals is a list of dictionaries. Get rid of the keys
          # so that it is a list of lists and sort each list in
          # the increasing order of base[l]^x
          for l, t in enumerate(transversals):
              transversals[l] = sorted(t.values(),
                                  key = lambda x: base_ordering[base[l]^x])
  
          orbits = H.basic_orbits
          h_stabs = H.basic_stabilizers
          g_stabs = self.basic_stabilizers
  
          indices = [x.order()//y.order() for x, y in zip(g_stabs, h_stabs)]
  
          # T^(l) should be a right transversal of H^(l) in G^(l) for
          # 1<=l<=len(base). While H^(l) is the trivial group, T^(l)
          # contains all the elements of G^(l) so we might just as well
          # start with l = len(h_stabs)-1
          T = g_stabs[len(h_stabs)]._elements
          t_len = len(T)
          l = len(h_stabs)-1
          while l > -1:
              T_next = []
              for u in transversals[l]:
                  if u == identity:
                      continue
                  b = base_ordering[base[l]^u]
                  for t in T:
                      p = t*u
                      if all([base_ordering[h^p] >= b for h in orbits[l]]):
                          T_next.append(p)
                      if t_len + len(T_next) == indices[l]:
                          break
                  if t_len + len(T_next) == indices[l]:
                      break
              T += T_next
              t_len += len(T_next)
              l -= 1
          return T
  
      def _coset_representative(self, g, H):
          """Return the representative of Hg from the transversal that
          would be computed by `self.coset_transversal(H)`.
  
          """
          if H.order() == 1:
              return g
          # The base of self must be an extension of H.base.
          if not(self.base[:len(H.base)] == H.base):
              self._schreier_sims(base=H.base)
          orbits = H.basic_orbits[:]
          h_transversals = [list(_.values()) for _ in H.basic_transversals]
          transversals = [list(_.values()) for _ in self.basic_transversals]
          base = self.base
          base_ordering = _base_ordering(base, self.degree)
          def step(l, x):
              gamma = sorted(orbits[l], key = lambda y: base_ordering[y^x])[0]
              i = [base[l]^h for h in h_transversals[l]].index(gamma)
              x = h_transversals[l][i]*x
              if l < len(orbits)-1:
                  for u in transversals[l]:
                      if base[l]^u == base[l]^x:
                          break
                  x = step(l+1, x*u**-1)*u
              return x
          return step(0, g)
  
      def coset_table(self, H):
          """Return the standardised (right) coset table of self in H as
          a list of lists.
          """
          # Maybe this should be made to return an instance of CosetTable
          # from fp_groups.py but the class would need to be changed first
          # to be compatible with PermutationGroups
  
          from itertools import chain, product
          if not H.is_subgroup(self):
              raise ValueError("The argument must be a subgroup")
          T = self.coset_transversal(H)
          n = len(T)
  
          A = list(chain.from_iterable((gen, gen**-1)
                      for gen in self.generators))
  
          table = []
          for i in range(n):
              row = [self._coset_representative(T[i]*x, H) for x in A]
              row = [T.index(r) for r in row]
              table.append(row)
  
  
          # standardize (this is the same as the algorithm used in fp_groups)
          # If CosetTable is made compatible with PermutationGroups, this
          # should be replaced by table.standardize()
          A = range(len(A))
          gamma = 1
          for alpha, a in product(range(n), A):
              beta = table[alpha][a]
              if beta >= gamma:
                  if beta > gamma:
                      for x in A:
                          z = table[gamma][x]
                          table[gamma][x] = table[beta][x]
                          table[beta][x] = z
                          for i in range(n):
                              if table[i][x] == beta:
                                  table[i][x] = gamma
                              elif table[i][x] == gamma:
                                  table[i][x] = beta
                  gamma += 1
                  if gamma == n-1:
                      return table
  
      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()``
  
          """
          return self.centralizer(self)
  
      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.
  
          Parameters
          ==========
  
          other
              a permutation group/list of permutations/single permutation
  
          Examples
          ========
  
          >>> from sympy.combinatorics.named_groups import (SymmetricGroup,
          ... CyclicGroup)
          >>> S = SymmetricGroup(6)
          >>> C = CyclicGroup(6)
          >>> H = S.centralizer(C)
          >>> H.is_subgroup(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 = _af_new(list(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]
                          im = g._array_form[base[l]]
                          im_rep = g._array_form[rep]
                          tr_el = transversals[rep_orb_index][base[l]]
                          # using the definition of transversal,
                          # base[l]^g = rep^(tr_el*g);
                          # if g belongs to the centralizer, then
                          # base[l]^g = (rep^g)^tr_el
                          return im == tr_el._array_form[im_rep]
                      tests[l] = test
  
              def prop(g):
                  return [rmul(g, gen) for gen in other.generators] == \
                         [rmul(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.is_subgroup(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)
  
          """
          ggens = G.generators
          hgens = H.generators
          commutators = []
          for ggen in ggens:
              for hgen in hgens:
                  commutator = rmul(hgen, ggen, ~hgen, ~ggen)
                  if commutator not in commutators:
                      commutators.append(commutator)
          res = self.normal_closure(commutators)
          return res
  
      def coset_factor(self, g, factor_index=False):
          """Return ``G``'s (self's) coset factorization of ``g``
  
          If ``g`` is an element of ``G`` then it can be written as the product
          of permutations drawn from the Schreier-Sims coset decomposition,
  
          The permutations returned in ``f`` are those for which
          the product gives ``g``: ``g = f[n]*...f[1]*f[0]`` where ``n = len(B)``
          and ``B = G.base``. f[i] is one of the permutations in
          ``self._basic_orbits[i]``.
  
          If factor_index==True,
          returns a tuple ``[b[0],..,b[n]]``, where ``b[i]``
          belongs to ``self._basic_orbits[i]``
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation, PermutationGroup
          >>> Permutation.print_cyclic = True
          >>> a = Permutation(0, 1, 3, 7, 6, 4)(2, 5)
          >>> b = Permutation(0, 1, 3, 2)(4, 5, 7, 6)
          >>> G = PermutationGroup([a, b])
  
          Define g:
  
          >>> g = Permutation(7)(1, 2, 4)(3, 6, 5)
  
          Confirm that it is an element of G:
  
          >>> G.contains(g)
          True
  
          Thus, it can be written as a product of factors (up to
          3) drawn from u. See below that a factor from u1 and u2
          and the Identity permutation have been used:
  
          >>> f = G.coset_factor(g)
          >>> f[2]*f[1]*f[0] == g
          True
          >>> f1 = G.coset_factor(g, True); f1
          [0, 4, 4]
          >>> tr = G.basic_transversals
          >>> f[0] == tr[0][f1[0]]
          True
  
          If g is not an element of G then [] is returned:
  
          >>> c = Permutation(5, 6, 7)
          >>> G.coset_factor(c)
          []
  
          see util._strip
          """
          if isinstance(g, (Cycle, Permutation)):
              g = g.list()
          if len(g) != self._degree:
              # this could either adjust the size or return [] immediately
              # but we don't choose between the two and just signal a possible
              # error
              raise ValueError('g should be the same size as permutations of G')
          I = list(range(self._degree))
          basic_orbits = self.basic_orbits
          transversals = self._transversals
          factors = []
          base = self.base
          h = g
          for i in range(len(base)):
              beta = h[base[i]]
              if beta == base[i]:
                  factors.append(beta)
                  continue
              if beta not in basic_orbits[i]:
                  return []
              u = transversals[i][beta]._array_form
              h = _af_rmul(_af_invert(u), h)
              factors.append(beta)
          if h != I:
              return []
          if factor_index:
              return factors
          tr = self.basic_transversals
          factors = [tr[i][factors[i]] for i in range(len(base))]
          return factors
  
      def generator_product(self, g):
          '''
          Return a list of strong generators `[s1, ..., sn]`
          s.t `g = sn*...*s1`.
  
          '''
          if g in self.strong_gens:
              return [g]
          f = self.coset_factor(g, True)
          product = []
          for i, j in enumerate(f):
              slp = self._transversal_slp[i][j]
              for s in slp:
                  product.append(self.strong_gens[s])
          return product
  
      def coset_rank(self, g):
          """rank using Schreier-Sims representation
  
          The coset rank of ``g`` is the ordering number in which
          it appears in the lexicographic listing according to the
          coset decomposition
  
          The ordering is the same as in G.generate(method='coset').
          If ``g`` does not belong to the group it returns None.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> a = Permutation(0, 1, 3, 7, 6, 4)(2, 5)
          >>> b = Permutation(0, 1, 3, 2)(4, 5, 7, 6)
          >>> G = PermutationGroup([a, b])
          >>> c = Permutation(7)(2, 4)(3, 5)
          >>> G.coset_rank(c)
          16
          >>> G.coset_unrank(16)
          (7)(2 4)(3 5)
  
          See Also
          ========
  
          coset_factor
  
          """
          factors = self.coset_factor(g, True)
          if not factors:
              return None
          rank = 0
          b = 1
          transversals = self._transversals
          base = self._base
          basic_orbits = self._basic_orbits
          for i in range(len(base)):
              k = factors[i]
              j = basic_orbits[i].index(k)
              rank += b*j
              b = b*len(transversals[i])
          return rank
  
      def coset_unrank(self, rank, af=False):
          """unrank using Schreier-Sims representation
  
          coset_unrank is the inverse operation of coset_rank
          if 0 <= rank < order; otherwise it returns None.
  
          """
          if rank < 0 or rank >= self.order():
              return None
          base = self._base
          transversals = self._transversals
          basic_orbits = self._basic_orbits
          m = len(base)
          v = [0]*m
          for i in range(m):
              rank, c = divmod(rank, len(transversals[i]))
              v[i] = basic_orbits[i][c]
          a = [transversals[i][v[i]]._array_form for i in range(m)]
          h = _af_rmuln(*a)
          if af:
              return h
          else:
              return _af_new(h)
  
      @property
      def degree(self):
          """Returns the size of the permutations in the group.
  
          The number of permutations comprising the group is given by
          ``len(group)``; the number of permutations that can be generated
          by the group is given by ``group.order()``.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> a = Permutation([1, 0, 2])
          >>> G = PermutationGroup([a])
          >>> G.degree
          3
          >>> len(G)
          1
          >>> G.order()
          2
          >>> list(G.generate())
          [(2), (2)(0 1)]
  
          See Also
          ========
  
          order
          """
          return self._degree
  
      @property
      def identity(self):
          '''
          Return the identity element of the permutation group.
  
          '''
          return _af_new(list(range(self.degree)))
  
      @property
      def elements(self):
          """Returns all the elements of the permutation group as a set
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation, PermutationGroup
          >>> p = PermutationGroup(Permutation(1, 3), Permutation(1, 2))
          >>> p.elements
          {(3), (2 3), (3)(1 2), (1 2 3), (1 3 2), (1 3)}
  
          """
          return set(self._elements)
  
      @property
      def _elements(self):
          """Returns all the elements of the permutation group as a list
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation, PermutationGroup
          >>> p = PermutationGroup(Permutation(1, 3), Permutation(1, 2))
          >>> p._elements
          [(3), (3)(1 2), (1 3), (2 3), (1 2 3), (1 3 2)]
  
          """
          return list(islice(self.generate(), None))
  
      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].is_subgroup(AlternatingGroup(4))
          True
          >>> S.derived_series()[2].is_subgroup(DihedralGroup(2))
          True
  
          See Also
          ========
  
          derived_subgroup
  
          """
          res = [self]
          current = self
          next = self.derived_subgroup()
          while not current.is_subgroup(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 import Permutation
          >>> Permutation.print_cyclic = True
          >>> 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 = [_af_invert(p) for p in gens]
          set_commutators = set()
          degree = self._degree
          rng = list(range(degree))
          for i in range(r):
              for j in range(r):
                  p1 = gens[i]
                  p2 = gens[j]
                  c = list(range(degree))
                  for k in rng:
                      c[p2[p1[k]]] = p1[p2[k]]
                  ct = tuple(c)
                  if not ct in set_commutators:
                      set_commutators.add(ct)
          cms = [_af_new(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
  
          Iteration is done with one of these methods::
  
            method='coset'  using the Schreier-Sims coset representation
            method='dimino' using the Dimino method
  
          If af = True it yields the array form of the permutations
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics import PermutationGroup
          >>> from sympy.combinatorics.polyhedron import tetrahedron
  
          The permutation group given in the tetrahedron object is also
          true groups:
  
          >>> G = tetrahedron.pgroup
          >>> G.is_group
          True
  
          Also the group generated by the permutations in the tetrahedron
          pgroup -- even the first two -- is a proper group:
  
          >>> H = PermutationGroup(G[0], G[1])
          >>> J = PermutationGroup(list(H.generate())); J
          PermutationGroup([
              (0 1)(2 3),
              (3),
              (1 2 3),
              (1 3 2),
              (0 3 1),
              (0 2 3),
              (0 3)(1 2),
              (0 1 3),
              (3)(0 2 1),
              (0 3 2),
              (3)(0 1 2),
              (0 2)(1 3)])
          >>> _.is_group
          True
          """
          if method == "coset":
              return self.generate_schreier_sims(af)
          elif method == "dimino":
              return self.generate_dimino(af)
          else:
              raise NotImplementedError('No generation defined for %s' % method)
  
      def generate_dimino(self, af=False):
          """Yield group elements using Dimino's algorithm
  
          If af == True it yields the array form of the permutations
  
          References
          ==========
  
          [1] The Implementation of Various Algorithms for Permutation Groups in
          the Computer Algebra System: AXIOM, N.J. Doye, M.Sc. Thesis
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> a = Permutation([0, 2, 1, 3])
          >>> b = Permutation([0, 2, 3, 1])
          >>> g = PermutationGroup([a, b])
          >>> list(g.generate_dimino(af=True))
          [[0, 1, 2, 3], [0, 2, 1, 3], [0, 2, 3, 1],
           [0, 1, 3, 2], [0, 3, 2, 1], [0, 3, 1, 2]]
  
          """
          idn = list(range(self.degree))
          order = 0
          element_list = [idn]
          set_element_list = {tuple(idn)}
          if af:
              yield idn
          else:
              yield _af_new(idn)
          gens = [p._array_form for p in self.generators]
  
          for i in range(len(gens)):
              # D elements of the subgroup G_i generated by gens[:i]
              D = element_list[:]
              N = [idn]
              while N:
                  A = N
                  N = []
                  for a in A:
                      for g in gens[:i + 1]:
                          ag = _af_rmul(a, g)
                          if tuple(ag) not in set_element_list:
                              # produce G_i*g
                              for d in D:
                                  order += 1
                                  ap = _af_rmul(d, ag)
                                  if af:
                                      yield ap
                                  else:
                                      p = _af_new(ap)
                                      yield p
                                  element_list.append(ap)
                                  set_element_list.add(tuple(ap))
                                  N.append(ap)
          self._order = len(element_list)
  
      def generate_schreier_sims(self, af=False):
          """Yield group elements using the Schreier-Sims representation
          in coset_rank order
  
          If ``af = True`` it yields the array form of the permutations
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> a = Permutation([0, 2, 1, 3])
          >>> b = Permutation([0, 2, 3, 1])
          >>> g = PermutationGroup([a, b])
          >>> list(g.generate_schreier_sims(af=True))
          [[0, 1, 2, 3], [0, 2, 1, 3], [0, 3, 2, 1],
           [0, 1, 3, 2], [0, 2, 3, 1], [0, 3, 1, 2]]
          """
  
          n = self._degree
          u = self.basic_transversals
          basic_orbits = self._basic_orbits
          if len(u) == 0:
              for x in self.generators:
                  if af:
                      yield x._array_form
                  else:
                      yield x
              return
          if len(u) == 1:
              for i in basic_orbits[0]:
                  if af:
                      yield u[0][i]._array_form
                  else:
                      yield u[0][i]
              return
  
          u = list(reversed(u))
          basic_orbits = basic_orbits[::-1]
          # stg stack of group elements
          stg = [list(range(n))]
          posmax = [len(x) for x in u]
          n1 = len(posmax) - 1
          pos = [0]*n1
          h = 0
          while 1:
              # backtrack when finished iterating over coset
              if pos[h] >= posmax[h]:
                  if h == 0:
                      return
                  pos[h] = 0
                  h -= 1
                  stg.pop()
                  continue
              p = _af_rmul(u[h][basic_orbits[h][pos[h]]]._array_form, stg[-1])
              pos[h] += 1
              stg.append(p)
              h += 1
              if h == n1:
                  if af:
                      for i in basic_orbits[-1]:
                          p = _af_rmul(u[-1][i]._array_form, stg[-1])
                          yield p
                  else:
                      for i in basic_orbits[-1]:
                          p = _af_rmul(u[-1][i]._array_form, stg[-1])
                          p1 = _af_new(p)
                          yield p1
                  stg.pop()
                  h -= 1
  
      @property
      def generators(self):
          """Returns the generators of the group.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> a = Permutation([0, 2, 1])
          >>> b = Permutation([1, 0, 2])
          >>> G = PermutationGroup([a, b])
          >>> G.generators
          [(1 2), (2)(0 1)]
  
          """
          return self._generators
  
      def contains(self, g, strict=True):
          """Test if permutation ``g`` belong to self, ``G``.
  
          If ``g`` is an element of ``G`` it can be written as a product
          of factors drawn from the cosets of ``G``'s stabilizers. To see
          if ``g`` is one of the actual generators defining the group use
          ``G.has(g)``.
  
          If ``strict`` is not ``True``, ``g`` will be resized, if necessary,
          to match the size of permutations in ``self``.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
  
          >>> a = Permutation(1, 2)
          >>> b = Permutation(2, 3, 1)
          >>> G = PermutationGroup(a, b, degree=5)
          >>> G.contains(G[0]) # trivial check
          True
          >>> elem = Permutation([[2, 3]], size=5)
          >>> G.contains(elem)
          True
          >>> G.contains(Permutation(4)(0, 1, 2, 3))
          False
  
          If strict is False, a permutation will be resized, if
          necessary:
  
          >>> H = PermutationGroup(Permutation(5))
          >>> H.contains(Permutation(3))
          False
          >>> H.contains(Permutation(3), strict=False)
          True
  
          To test if a given permutation is present in the group:
  
          >>> elem in G.generators
          False
          >>> G.has(elem)
          False
  
          See Also
          ========
  
          coset_factor, has, in
  
          """
          if not isinstance(g, Permutation):
              return False
          if g.size != self.degree:
              if strict:
                  return False
              g = Permutation(g, size=self.degree)
          if g in self.generators:
              return True
          return bool(self.coset_factor(g.array_form, True))
  
      @property
      def is_abelian(self):
          """Test if the group is Abelian.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> a = Permutation([0, 2, 1])
          >>> b = Permutation([1, 0, 2])
          >>> G = PermutationGroup([a, b])
          >>> G.is_abelian
          False
          >>> a = Permutation([0, 2, 1])
          >>> G = PermutationGroup([a])
          >>> G.is_abelian
          True
  
          """
          if self._is_abelian is not None:
              return self._is_abelian
  
          self._is_abelian = True
          gens = [p._array_form for p in self.generators]
          for x in gens:
              for y in gens:
                  if y <= x:
                      continue
                  if not _af_commutes_with(x, y):
                      self._is_abelian = False
                      return False
          return True
  
      def is_alt_sym(self, eps=0.05, _random_prec=None):
          r"""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.
  
          Notes
          =====
  
          The algorithm itself uses some nontrivial results from group theory and
          number theory:
          1) If a transitive group ``G`` of degree ``n`` contains an element
          with a cycle of length ``n/2 < p < n-2`` for ``p`` a prime, ``G`` is the
          symmetric or alternating group ([1], pp. 81-82)
          2) The proportion of elements in the symmetric/alternating group having
          the property described in 1) is approximately `\log(2)/\log(n)`
          ([1], p.82; [2], pp. 226-227).
          The helper function ``_check_cycles_alt_sym`` is used to
          go over the cycles in a permutation and look for ones satisfying 1).
  
          Examples
          ========
  
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> from sympy.combinatorics.named_groups import DihedralGroup
          >>> D = DihedralGroup(10)
          >>> D.is_alt_sym()
          False
  
          See Also
          ========
  
          _check_cycles_alt_sym
  
          """
          if _random_prec is None:
              n = self.degree
              if n < 8:
                  return False
              if not self.is_transitive():
                  return False
              if n < 17:
                  c_n = 0.34
              else:
                  c_n = 0.57
              d_n = (c_n*log(2))/log(n)
              N_eps = int(-log(eps)/d_n)
              for i in range(N_eps):
                  perm = self.random_pr()
                  if _check_cycles_alt_sym(perm):
                      return True
              return False
          else:
              for i in range(_random_prec['N_eps']):
                  perm = _random_prec[i]
                  if _check_cycles_alt_sym(perm):
                      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 = _af_new(list(range(degree)))
              if all(g == identity for g in 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, strict=True):
          """Test if ``G=self`` is a normal subgroup of ``gr``.
  
          G is normal in gr if
          for each g2 in G, g1 in gr, ``g = g1*g2*g1**-1`` belongs to G
          It is sufficient to check this for each g1 in gr.generators and
          g2 in G.generators.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> a = Permutation([1, 2, 0])
          >>> b = Permutation([1, 0, 2])
          >>> G = PermutationGroup([a, b])
          >>> G1 = PermutationGroup([a, Permutation([2, 0, 1])])
          >>> G1.is_normal(G)
          True
  
          """
          d_self = self.degree
          d_gr = gr.degree
          new_self = self.copy()
          if not strict and d_self != d_gr:
              if d_self < d_gr:
                  new_self = PermGroup(new_self.generators + [Permutation(d_gr - 1)])
              else:
                  gr = PermGroup(gr.generators + [Permutation(d_self - 1)])
          gens2 = [p._array_form for p in new_self.generators]
          gens1 = [p._array_form for p in gr.generators]
          for g1 in gens1:
              for g2 in gens2:
                  p = _af_rmuln(g1, g2, _af_invert(g1))
                  if not new_self.coset_factor(p, True):
                      return False
          return True
  
      def is_primitive(self, randomized=True):
          """Test if a group is primitive.
  
          A permutation group ``G`` acting on a set ``S`` is called primitive if
          ``S`` contains no nontrivial block under the action of ``G``
          (a block is nontrivial if its cardinality is more than ``1``).
  
          Notes
          =====
  
          The algorithm is described in [1], p.83, and uses the function
          minimal_block to search for blocks of the form `\{0, k\}` for ``k``
          ranging over representatives for the orbits of `G_0`, the stabilizer of
          ``0``. This algorithm has complexity `O(n^2)` where ``n`` is the degree
          of the group, and will perform badly if `G_0` is small.
  
          There are two implementations offered: one finds `G_0`
          deterministically using the function ``stabilizer``, and the other
          (default) produces random elements of `G_0` using ``random_stab``,
          hoping that they generate a subgroup of `G_0` with not too many more
          orbits than `G_0` (this is suggested in [1], p.83). Behavior is changed
          by the ``randomized`` flag.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> from sympy.combinatorics.named_groups import DihedralGroup
          >>> D = DihedralGroup(10)
          >>> D.is_primitive()
          False
  
          See Also
          ========
  
          minimal_block, random_stab
  
          """
          if self._is_primitive is not None:
              return self._is_primitive
          n = self.degree
          if randomized:
              random_stab_gens = []
              v = self.schreier_vector(0)
              for i in range(len(self)):
                  random_stab_gens.append(self.random_stab(0, v))
              stab = PermutationGroup(random_stab_gens)
          else:
              stab = self.stabilizer(0)
          orbits = stab.orbits()
          for orb in orbits:
              x = orb.pop()
              if x != 0 and self.minimal_block([0, x]) != [0]*n:
                  self._is_primitive = False
                  return False
          self._is_primitive = True
          return True
  
      @property
      def is_solvable(self):
          """Test if the group is solvable.
  
          ``G`` is solvable if its derived series terminates with the trivial
          group ([1], p.29).
  
          Examples
          ========
  
          >>> from sympy.combinatorics.named_groups import SymmetricGroup
          >>> S = SymmetricGroup(3)
          >>> S.is_solvable
          True
  
          See Also
          ========
  
          is_nilpotent, derived_series
  
          """
          if self._is_solvable is None:
              ds = self.derived_series()
              terminator = ds[len(ds) - 1]
              gens = terminator.generators
              degree = self.degree
              identity = _af_new(list(range(degree)))
              if all(g == identity for g in gens):
                  self._is_solvable = True
                  return True
              else:
                  self._is_solvable = False
                  return False
          else:
              return self._is_solvable
  
      def is_subgroup(self, G, strict=True):
          """Return ``True`` if all elements of ``self`` belong to ``G``.
  
          If ``strict`` is ``False`` then if ``self``'s degree is smaller
          than ``G``'s, the elements will be resized to have the same degree.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation, PermutationGroup
          >>> from sympy.combinatorics.named_groups import (SymmetricGroup,
          ...    CyclicGroup)
  
          Testing is strict by default: the degree of each group must be the
          same:
  
          >>> p = Permutation(0, 1, 2, 3, 4, 5)
          >>> G1 = PermutationGroup([Permutation(0, 1, 2), Permutation(0, 1)])
          >>> G2 = PermutationGroup([Permutation(0, 2), Permutation(0, 1, 2)])
          >>> G3 = PermutationGroup([p, p**2])
          >>> assert G1.order() == G2.order() == G3.order() == 6
          >>> G1.is_subgroup(G2)
          True
          >>> G1.is_subgroup(G3)
          False
          >>> G3.is_subgroup(PermutationGroup(G3[1]))
          False
          >>> G3.is_subgroup(PermutationGroup(G3[0]))
          True
  
          To ignore the size, set ``strict`` to ``False``:
  
          >>> S3 = SymmetricGroup(3)
          >>> S5 = SymmetricGroup(5)
          >>> S3.is_subgroup(S5, strict=False)
          True
          >>> C7 = CyclicGroup(7)
          >>> G = S5*C7
          >>> S5.is_subgroup(G, False)
          True
          >>> C7.is_subgroup(G, 0)
          False
  
          """
          if not isinstance(G, PermutationGroup):
              return False
          if self == G:
              return True
          if G.order() % self.order() != 0:
              return False
          if self.degree == G.degree or \
                  (self.degree < G.degree and not strict):
              gens = self.generators
          else:
              return False
          return all(G.contains(g, strict=strict) for g in gens)
  
      def is_transitive(self, strict=True):
          """Test if the group is transitive.
  
          A group is transitive if it has a single orbit.
  
          If ``strict`` is ``False`` the group is transitive if it has
          a single orbit of length different from 1.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.permutations import Permutation
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> a = Permutation([0, 2, 1, 3])
          >>> b = Permutation([2, 0, 1, 3])
          >>> G1 = PermutationGroup([a, b])
          >>> G1.is_transitive()
          False
          >>> G1.is_transitive(strict=False)
          True
          >>> c = Permutation([2, 3, 0, 1])
          >>> G2 = PermutationGroup([a, c])
          >>> G2.is_transitive()
          True
          >>> d = Permutation([1, 0, 2, 3])
          >>> e = Permutation([0, 1, 3, 2])
          >>> G3 = PermutationGroup([d, e])
          >>> G3.is_transitive() or G3.is_transitive(strict=False)
          False
  
          """
          if self._is_transitive:  # strict or not, if True then True
              return self._is_transitive
          if strict:
              if self._is_transitive is not None:  # we only store strict=True
                  return self._is_transitive
  
              ans = len(self.orbit(0)) == self.degree
              self._is_transitive = ans
              return ans
  
          got_orb = False
          for x in self.orbits():
              if len(x) > 1:
                  if got_orb:
                      return False
                  got_orb = True
          return got_orb
  
      @property
      def is_trivial(self):
          """Test if the group is the trivial group.
  
          This is true if the group contains only the identity permutation.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> G = PermutationGroup([Permutation([0, 1, 2])])
          >>> G.is_trivial
          True
  
          """
          if self._is_trivial is None:
              self._is_trivial = len(self) == 1 and self[0].is_Identity
          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].is_subgroup(DihedralGroup(2))
          True
  
          See Also
          ========
  
          commutator, derived_series
  
          """
          res = [self]
          current = self
          next = self.commutator(self, current)
          while not current.is_subgroup(next):
              res.append(next)
              current = next
              next = self.commutator(self, current)
          return res
  
      @property
      def max_div(self):
          """Maximum proper divisor of the degree of a permutation group.
  
          Notes
          =====
  
          Obviously, this is the degree divided by its minimal proper divisor
          (larger than ``1``, if one exists). As it is guaranteed to be prime,
          the ``sieve`` from ``sympy.ntheory`` is used.
          This function is also used as an optimization tool for the functions
          ``minimal_block`` and ``_union_find_merge``.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> G = PermutationGroup([Permutation([0, 2, 1, 3])])
          >>> G.max_div
          2
  
          See Also
          ========
  
          minimal_block, _union_find_merge
  
          """
          if self._max_div is not None:
              return self._max_div
          n = self.degree
          if n == 1:
              return 1
          for x in sieve:
              if n % x == 0:
                  d = n//x
                  self._max_div = d
                  return d
  
      def minimal_block(self, points):
          r"""For a transitive group, finds the block system generated by
          ``points``.
  
          If a group ``G`` acts on a set ``S``, a nonempty subset ``B`` of ``S``
          is called a block under the action of ``G`` if for all ``g`` in ``G``
          we have ``gB = B`` (``g`` fixes ``B``) or ``gB`` and ``B`` have no
          common points (``g`` moves ``B`` entirely). ([1], p.23; [6]).
  
          The distinct translates ``gB`` of a block ``B`` for ``g`` in ``G``
          partition the set ``S`` and this set of translates is known as a block
          system. Moreover, we obviously have that all blocks in the partition
          have the same size, hence the block size divides ``|S|`` ([1], p.23).
          A ``G``-congruence is an equivalence relation ``~`` on the set ``S``
          such that ``a ~ b`` implies ``g(a) ~ g(b)`` for all ``g`` in ``G``.
          For a transitive group, the equivalence classes of a ``G``-congruence
          and the blocks of a block system are the same thing ([1], p.23).
  
          The algorithm below checks the group for transitivity, and then finds
          the ``G``-congruence generated by the pairs ``(p_0, p_1), (p_0, p_2),
          ..., (p_0,p_{k-1})`` which is the same as finding the maximal block
          system (i.e., the one with minimum block size) such that
          ``p_0, ..., p_{k-1}`` are in the same block ([1], p.83).
  
          It is an implementation of Atkinson's algorithm, as suggested in [1],
          and manipulates an equivalence relation on the set ``S`` using a
          union-find data structure. The running time is just above
          `O(|points||S|)`. ([1], pp. 83-87; [7]).
  
          Examples
          ========
  
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> from sympy.combinatorics.named_groups import DihedralGroup
          >>> D = DihedralGroup(10)
          >>> D.minimal_block([0, 5])
          [0, 6, 2, 8, 4, 0, 6, 2, 8, 4]
          >>> D.minimal_block([0, 1])
          [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
  
          See Also
          ========
  
          _union_find_rep, _union_find_merge, is_transitive, is_primitive
  
          """
          if not self.is_transitive():
              return False
          n = self.degree
          gens = self.generators
          # initialize the list of equivalence class representatives
          parents = list(range(n))
          ranks = [1]*n
          not_rep = []
          k = len(points)
          # the block size must divide the degree of the group
          if k > self.max_div:
              return [0]*n
          for i in range(k - 1):
              parents[points[i + 1]] = points[0]
              not_rep.append(points[i + 1])
          ranks[points[0]] = k
          i = 0
          len_not_rep = k - 1
          while i < len_not_rep:
              gamma = not_rep[i]
              i += 1
              for gen in gens:
                  # find has side effects: performs path compression on the list
                  # of representatives
                  delta = self._union_find_rep(gamma, parents)
                  # union has side effects: performs union by rank on the list
                  # of representatives
                  temp = self._union_find_merge(gen(gamma), gen(delta), ranks,
                                                parents, not_rep)
                  if temp == -1:
                      return [0]*n
                  len_not_rep += temp
          for i in range(n):
              # force path compression to get the final state of the equivalence
              # relation
              self._union_find_rep(i, parents)
          return parents
  
      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.named_groups import (SymmetricGroup,
          ... CyclicGroup, AlternatingGroup)
          >>> S = SymmetricGroup(5)
          >>> C = CyclicGroup(5)
          >>> G = S.normal_closure(C)
          >>> G.order()
          60
          >>> G.is_subgroup(AlternatingGroup(5))
          True
  
          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 = _af_new(list(range(degree)))
  
              if all(g == identity for g 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)
  
              self._random_pr_init(r=10, n=20)
  
              _loop = True
              while _loop:
                  Z._random_pr_init(r=10, n=10)
                  for i in range(k):
                      g = self.random_pr()
                      h = Z.random_pr()
                      conj = 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)
                  _loop = False
                  for g in self.generators:
                      for h in Z.generators:
                          conj = h^g
                          res = _strip(conj, base, basic_orbits,
                                       basic_transversals)
                          if res[0] != identity or res[1] != len(base) + 1:
                              _loop = True
                              break
                      if _loop:
                          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"""Compute the orbit of alpha `\{g(\alpha) | g \in G\}` as a set.
  
          The time complexity of the algorithm used here is `O(|Orb|*r)` where
          `|Orb|` is the size of the orbit and ``r`` is the number of generators of
          the group. For a more detailed analysis, see [1], p.78, [2], pp. 19-21.
          Here alpha can be a single point, or a list of points.
  
          If alpha is a single point, the ordinary orbit is computed.
          if alpha is a list of points, there are three available options:
  
          'union' - computes the union of the orbits of the points in the list
          'tuples' - computes the orbit of the list interpreted as an ordered
          tuple under the group action ( i.e., g((1,2,3)) = (g(1), g(2), g(3)) )
          'sets' - computes the orbit of the list interpreted as a sets
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> a = Permutation([1, 2, 0, 4, 5, 6, 3])
          >>> G = PermutationGroup([a])
          >>> G.orbit(0)
          {0, 1, 2}
          >>> G.orbit([0, 4], 'union')
          {0, 1, 2, 3, 4, 5, 6}
  
          See Also
          ========
  
          orbit_transversal
  
          """
          return _orbit(self.degree, self.generators, alpha, action)
  
      def orbit_rep(self, alpha, beta, schreier_vector=None):
          """Return a group element which sends ``alpha`` to ``beta``.
  
          If ``beta`` is not in the orbit of ``alpha``, the function returns
          ``False``. This implementation makes use of the schreier vector.
          For a proof of correctness, see [1], p.80
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> from sympy.combinatorics.named_groups import AlternatingGroup
          >>> G = AlternatingGroup(5)
          >>> G.orbit_rep(0, 4)
          (0 4 1 2 3)
  
          See Also
          ========
  
          schreier_vector
  
          """
          if schreier_vector is None:
              schreier_vector = self.schreier_vector(alpha)
          if schreier_vector[beta] is None:
              return False
          k = schreier_vector[beta]
          gens = [x._array_form for x in self.generators]
          a = []
          while k != -1:
              a.append(gens[k])
              beta = gens[k].index(beta) # beta = (~gens[k])(beta)
              k = schreier_vector[beta]
          if a:
              return _af_new(_af_rmuln(*a))
          else:
              return _af_new(list(range(self._degree)))
  
      def orbit_transversal(self, alpha, pairs=False):
          r"""Computes a transversal for the orbit of ``alpha`` as a set.
  
          For a permutation group `G`, a transversal for the orbit
          `Orb = \{g(\alpha) | g \in G\}` is a set
          `\{g_\beta | g_\beta(\alpha) = \beta\}` for `\beta \in Orb`.
          Note that there may be more than one possible transversal.
          If ``pairs`` is set to ``True``, it returns the list of pairs
          `(\beta, g_\beta)`. For a proof of correctness, see [1], p.79
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> from sympy.combinatorics.named_groups import DihedralGroup
          >>> G = DihedralGroup(6)
          >>> G.orbit_transversal(0)
          [(5), (0 1 2 3 4 5), (0 5)(1 4)(2 3), (0 2 4)(1 3 5), (5)(0 4)(1 3), (0 3)(1 4)(2 5)]
  
          See Also
          ========
  
          orbit
  
          """
          return _orbit_transversal(self._degree, self.generators, alpha, pairs)
  
      def orbits(self, rep=False):
          """Return the orbits of ``self``, ordered according to lowest element
          in each orbit.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.permutations import Permutation
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> a = Permutation(1, 5)(2, 3)(4, 0, 6)
          >>> b = Permutation(1, 5)(3, 4)(2, 6, 0)
          >>> G = PermutationGroup([a, b])
          >>> G.orbits()
          [{0, 2, 3, 4, 6}, {1, 5}]
          """
          return _orbits(self._degree, self._generators)
  
      def order(self):
          """Return the order of the group: the number of permutations that
          can be generated from elements of the group.
  
          The number of permutations comprising the group is given by
          ``len(group)``; the length of each permutation in the group is
          given by ``group.size``.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.permutations import Permutation
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
  
          >>> a = Permutation([1, 0, 2])
          >>> G = PermutationGroup([a])
          >>> G.degree
          3
          >>> len(G)
          1
          >>> G.order()
          2
          >>> list(G.generate())
          [(2), (2)(0 1)]
  
          >>> a = Permutation([0, 2, 1])
          >>> b = Permutation([1, 0, 2])
          >>> G = PermutationGroup([a, b])
          >>> G.order()
          6
  
          See Also
          ========
          degree
  
          """
          if self._order != None:
              return self._order
          if self._is_sym:
              n = self._degree
              self._order = factorial(n)
              return self._order
          if self._is_alt:
              n = self._degree
              self._order = factorial(n)/2
              return self._order
  
          basic_transversals = self.basic_transversals
          m = 1
          for x in basic_transversals:
              m *= len(x)
          self._order = m
          return m
  
      def pointwise_stabilizer(self, points, incremental=True):
          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.is_subgroup(S.stabilizer(2).stabilizer(3).stabilizer(5))
          True
  
          See Also
          ========
  
          stabilizer, schreier_sims_incremental
  
          Notes
          =====
  
          When incremental == True,
          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.
  
          """
          if incremental:
              base, strong_gens = self.schreier_sims_incremental(base=points)
              stab_gens = []
              degree = self.degree
              for gen in strong_gens:
                  if [gen(point) for point in points] == points:
                      stab_gens.append(gen)
              if not stab_gens:
                  stab_gens = _af_new(list(range(degree)))
              return PermutationGroup(stab_gens)
          else:
              gens = self._generators
              degree = self.degree
              for x in points:
                  gens = _stabilizer(degree, gens, x)
          return PermutationGroup(gens)
  
      def make_perm(self, n, seed=None):
          """
          Multiply ``n`` randomly selected permutations from
          pgroup together, starting with the identity
          permutation. If ``n`` is a list of integers, those
          integers will be used to select the permutations and they
          will be applied in L to R order: make_perm((A, B, C)) will
          give CBA(I) where I is the identity permutation.
  
          ``seed`` is used to set the seed for the random selection
          of permutations from pgroup. If this is a list of integers,
          the corresponding permutations from pgroup will be selected
          in the order give. This is mainly used for testing purposes.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> a, b = [Permutation([1, 0, 3, 2]), Permutation([1, 3, 0, 2])]
          >>> G = PermutationGroup([a, b])
          >>> G.make_perm(1, [0])
          (0 1)(2 3)
          >>> G.make_perm(3, [0, 1, 0])
          (0 2 3 1)
          >>> G.make_perm([0, 1, 0])
          (0 2 3 1)
  
          See Also
          ========
  
          random
          """
          if is_sequence(n):
              if seed is not None:
                  raise ValueError('If n is a sequence, seed should be None')
              n, seed = len(n), n
          else:
              try:
                  n = int(n)
              except TypeError:
                  raise ValueError('n must be an integer or a sequence.')
          randrange = _randrange(seed)
  
          # start with the identity permutation
          result = Permutation(list(range(self.degree)))
          m = len(self)
          for i in range(n):
              p = self[randrange(m)]
              result = rmul(result, p)
          return result
  
      def random(self, af=False):
          """Return a random group element
          """
          rank = randrange(self.order())
          return self.coset_unrank(rank, af)
  
      def random_pr(self, gen_count=11, iterations=50, _random_prec=None):
          """Return a random group element using product replacement.
  
          For the details of the product replacement algorithm, see
          ``_random_pr_init`` In ``random_pr`` the actual 'product replacement'
          is performed. Notice that if the attribute ``_random_gens``
          is empty, it needs to be initialized by ``_random_pr_init``.
  
          See Also
          ========
  
          _random_pr_init
  
          """
          if self._random_gens == []:
              self._random_pr_init(gen_count, iterations)
          random_gens = self._random_gens
          r = len(random_gens) - 1
  
          # handle randomized input for testing purposes
          if _random_prec is None:
              s = randrange(r)
              t = randrange(r - 1)
              if t == s:
                  t = r - 1
              x = choice([1, 2])
              e = choice([-1, 1])
          else:
              s = _random_prec['s']
              t = _random_prec['t']
              if t == s:
                  t = r - 1
              x = _random_prec['x']
              e = _random_prec['e']
  
          if x == 1:
              random_gens[s] = _af_rmul(random_gens[s], _af_pow(random_gens[t], e))
              random_gens[r] = _af_rmul(random_gens[r], random_gens[s])
          else:
              random_gens[s] = _af_rmul(_af_pow(random_gens[t], e), random_gens[s])
              random_gens[r] = _af_rmul(random_gens[s], random_gens[r])
          return _af_new(random_gens[r])
  
      def random_stab(self, alpha, schreier_vector=None, _random_prec=None):
          """Random element from the stabilizer of ``alpha``.
  
          The schreier vector for ``alpha`` is an optional argument used
          for speeding up repeated calls. The algorithm is described in [1], p.81
  
          See Also
          ========
  
          random_pr, orbit_rep
  
          """
          if schreier_vector is None:
              schreier_vector = self.schreier_vector(alpha)
          if _random_prec is None:
              rand = self.random_pr()
          else:
              rand = _random_prec['rand']
          beta = rand(alpha)
          h = self.orbit_rep(alpha, beta, schreier_vector)
          return rmul(~h, rand)
  
      def schreier_sims(self):
          """Schreier-Sims algorithm.
  
          It computes the generators of the chain of stabilizers
          `G > G_{b_1} > .. > G_{b1,..,b_r} > 1`
          in which `G_{b_1,..,b_i}` stabilizes `b_1,..,b_i`,
          and the corresponding ``s`` cosets.
          An element of the group can be written as the product
          `h_1*..*h_s`.
  
          We use the incremental Schreier-Sims algorithm.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.permutations import Permutation
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> a = Permutation([0, 2, 1])
          >>> b = Permutation([1, 0, 2])
          >>> G = PermutationGroup([a, b])
          >>> G.schreier_sims()
          >>> G.basic_transversals
          [{0: (2)(0 1), 1: (2), 2: (1 2)},
           {0: (2), 2: (0 2)}]
          """
          if self._transversals:
              return
          self._schreier_sims()
          return
  
      def _schreier_sims(self, base=None):
          base, strong_gens = self.schreier_sims_incremental(base=base)
          self._base = base
          self._strong_gens = strong_gens
          if not base:
              self._transversals = []
              self._basic_orbits = []
              return
  
          strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
          basic_orbits, transversals, slps = _orbits_transversals_from_bsgs(base,\
                  strong_gens_distr, slp=True)
  
          # rewrite the indices stored in slps in terms of strong_gens
          for i, slp in enumerate(slps):
              gens = strong_gens_distr[i]
              for k in slp:
                  slp[k] = [strong_gens.index(gens[s]) for s in slp[k]]
  
          self._transversals = transversals
          self._basic_orbits = [sorted(x) for x in basic_orbits]
          self._transversal_slp = slps
  
      def schreier_sims_incremental(self, base=None, gens=None):
          """Extend a sequence of points and generating set to a base and strong
          generating set.
  
          Parameters
          ==========
  
          base
              The sequence of points to be extended to a base. Optional
              parameter with default value ``[]``.
          gens
              The generating set to be extended to a strong generating set
              relative to the base obtained. Optional parameter with default
              value ``self.generators``.
  
          Returns
          =======
  
          (base, strong_gens)
              ``base`` is the base obtained, and ``strong_gens`` is the strong
              generating set relative to it. The original parameters ``base``,
              ``gens`` remain unchanged.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.named_groups import AlternatingGroup
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> from sympy.combinatorics.testutil import _verify_bsgs
          >>> A = AlternatingGroup(7)
          >>> base = [2, 3]
          >>> seq = [2, 3]
          >>> base, strong_gens = A.schreier_sims_incremental(base=seq)
          >>> _verify_bsgs(A, base, strong_gens)
          True
          >>> base[:2]
          [2, 3]
  
          Notes
          =====
  
          This version of the Schreier-Sims algorithm runs in polynomial time.
          There are certain assumptions in the implementation - if the trivial
          group is provided, ``base`` and ``gens`` are returned immediately,
          as any sequence of points is a base for the trivial group. If the
          identity is present in the generators ``gens``, it is removed as
          it is a redundant generator.
          The implementation is described in [1], pp. 90-93.
  
          See Also
          ========
  
          schreier_sims, schreier_sims_random
  
          """
          if base is None:
              base = []
          if gens is None:
              gens = self.generators[:]
          degree = self.degree
          id_af = list(range(degree))
          # handle the trivial group
          if len(gens) == 1 and gens[0].is_Identity:
              self._strong_gens_slp = {gens[0]: [gens[0]]}
              return base, gens
          # prevent side effects
          _base, _gens = base[:], gens[:]
          # remove the identity as a generator
          _gens = [x for x in _gens if not x.is_Identity]
          # make sure no generator fixes all base points
          for gen in _gens:
              if all(x == gen._array_form[x] for x in _base):
                  for new in id_af:
                      if gen._array_form[new] != new:
                          break
                  else:
                      assert None  # can this ever happen?
                  _base.append(new)
          # distribute generators according to basic stabilizers
          strong_gens_distr = _distribute_gens_by_base(_base, _gens)
          strong_gens_slp = {}
          # initialize the basic stabilizers, basic orbits and basic transversals
          orbs = {}
          transversals = {}
          slps = {}
          base_len = len(_base)
          for i in range(base_len):
              transversals[i], slps[i] = _orbit_transversal(degree, strong_gens_distr[i],
                  _base[i], pairs=True, af=True, slp=True)
              transversals[i] = dict(transversals[i])
              orbs[i] = list(transversals[i].keys())
          # main loop: amend the stabilizer chain until we have generators
          # for all stabilizers
          i = base_len - 1
          while i >= 0:
              # this flag is used to continue with the main loop from inside
              # a nested loop
              continue_i = False
              # test the generators for being a strong generating set
              db = {}
              for beta, u_beta in list(transversals[i].items()):
                  for j, gen in enumerate(strong_gens_distr[i]):
                      gb = gen._array_form[beta]
                      u1 = transversals[i][gb]
                      g1 = _af_rmul(gen._array_form, u_beta)
                      slp = [(i, g) for g in slps[i][beta]]
                      slp = [(i, j)] + slp
                      if g1 != u1:
                          # test if the schreier generator is in the i+1-th
                          # would-be basic stabilizer
                          y = True
                          try:
                              u1_inv = db[gb]
                          except KeyError:
                              u1_inv = db[gb] = _af_invert(u1)
                          schreier_gen = _af_rmul(u1_inv, g1)
                          u1_inv_slp = slps[i][gb][:]
                          u1_inv_slp.reverse()
                          u1_inv_slp = [(i, (g,)) for g in u1_inv_slp]
                          slp = u1_inv_slp + slp
                          h, j, slp = _strip_af(schreier_gen, _base, orbs, transversals, i, slp=slp, slps=slps)
                          if j <= base_len:
                              # new strong generator h at level j
                              y = False
                          elif h:
                              # h fixes all base points
                              y = False
                              moved = 0
                              while h[moved] == moved:
                                  moved += 1
                              _base.append(moved)
                              base_len += 1
                              strong_gens_distr.append([])
                          if y is False:
                              # if a new strong generator is found, update the
                              # data structures and start over
                              h = _af_new(h)
                              strong_gens_slp[h] = slp
                              for l in range(i + 1, j):
                                  strong_gens_distr[l].append(h)
                                  transversals[l], slps[l] =\
                                  _orbit_transversal(degree, strong_gens_distr[l],
                                      _base[l], pairs=True, af=True, slp=True)
                                  transversals[l] = dict(transversals[l])
                                  orbs[l] = list(transversals[l].keys())
                              i = j - 1
                              # continue main loop using the flag
                              continue_i = True
                      if continue_i is True:
                          break
                  if continue_i is True:
                      break
              if continue_i is True:
                  continue
              i -= 1
          # build the strong generating set
          strong_gens = list(uniq(i for gens in strong_gens_distr for i in gens))
  
          # rewrite the indices of strong_gens_slp in terms of the elements
          # of strong_gens
          for k in strong_gens_slp:
              slp = strong_gens_slp[k]
              for i in range(len(slp)):
                  s = slp[i]
                  if isinstance(s[1], tuple):
                      slp[i] = strong_gens_distr[s[0]][s[1][0]]**-1
                  else:
                      slp[i] = strong_gens_distr[s[0]][s[1]]
              strong_gens_slp[k] = slp
          # add the original generators
          for g in _gens:
              strong_gens_slp[g] = [strong_gens.index(g)]
          self._strong_gens_slp = strong_gens_slp
          return _base, strong_gens
  
      def schreier_sims_random(self, base=None, gens=None, consec_succ=10,
                               _random_prec=None):
          r"""Randomized Schreier-Sims algorithm.
  
          The randomized Schreier-Sims algorithm takes the sequence ``base``
          and the generating set ``gens``, and extends ``base`` to a base, and
          ``gens`` to a strong generating set relative to that base with
          probability of a wrong answer at most `2^{-consec\_succ}`,
          provided the random generators are sufficiently random.
  
          Parameters
          ==========
  
          base
              The sequence to be extended to a base.
          gens
              The generating set to be extended to a strong generating set.
          consec_succ
              The parameter defining the probability of a wrong answer.
          _random_prec
              An internal parameter used for testing purposes.
  
          Returns
          =======
  
          (base, strong_gens)
              ``base`` is the base and ``strong_gens`` is the strong generating
              set relative to it.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> 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)
          >>> _verify_bsgs(S, base, strong_gens) #doctest: +SKIP
          True
  
          Notes
          =====
  
          The algorithm is described in detail in [1], pp. 97-98. It extends
          the orbits ``orbs`` and the permutation groups ``stabs`` to
          basic orbits and basic stabilizers for the base and strong generating
          set produced in the end.
          The idea of the extension process
          is to "sift" random group elements through the stabilizer chain
          and amend the stabilizers/orbits along the way when a sift
          is not successful.
          The helper function ``_strip`` is used to attempt
          to decompose a random group element according to the current
          state of the stabilizer chain and report whether the element was
          fully decomposed (successful sift) or not (unsuccessful sift). In
          the latter case, the level at which the sift failed is reported and
          used to amend ``stabs``, ``base``, ``gens`` and ``orbs`` accordingly.
          The halting condition is for ``consec_succ`` consecutive successful
          sifts to pass. This makes sure that the current ``base`` and ``gens``
          form a BSGS with probability at least `1 - 1/\text{consec\_succ}`.
  
          See Also
          ========
  
          schreier_sims
  
          """
          if base is None:
              base = []
          if gens is None:
              gens = self.generators
          base_len = len(base)
          n = self.degree
          # make sure no generator fixes all base points
          for gen in gens:
              if all(gen(x) == x for x in base):
                  new = 0
                  while gen._array_form[new] == new:
                      new += 1
                  base.append(new)
                  base_len += 1
          # distribute generators according to basic stabilizers
          strong_gens_distr = _distribute_gens_by_base(base, gens)
          # initialize the basic stabilizers, basic transversals and basic orbits
          transversals = {}
          orbs = {}
          for i in range(base_len):
              transversals[i] = dict(_orbit_transversal(n, strong_gens_distr[i],
                  base[i], pairs=True))
              orbs[i] = list(transversals[i].keys())
          # initialize the number of consecutive elements sifted
          c = 0
          # start sifting random elements while the number of consecutive sifts
          # is less than consec_succ
          while c < consec_succ:
              if _random_prec is None:
                  g = self.random_pr()
              else:
                  g = _random_prec['g'].pop()
              h, j = _strip(g, base, orbs, transversals)
              y = True
              # determine whether a new base point is needed
              if j <= base_len:
                  y = False
              elif not h.is_Identity:
                  y = False
                  moved = 0
                  while h(moved) == moved:
                      moved += 1
                  base.append(moved)
                  base_len += 1
                  strong_gens_distr.append([])
              # if the element doesn't sift, amend the strong generators and
              # associated stabilizers and orbits
              if y is False:
                  for l in range(1, j):
                      strong_gens_distr[l].append(h)
                      transversals[l] = dict(_orbit_transversal(n,
                          strong_gens_distr[l], base[l], pairs=True))
                      orbs[l] = list(transversals[l].keys())
                  c = 0
              else:
                  c += 1
          # build the strong generating set
          strong_gens = strong_gens_distr[0][:]
          for gen in strong_gens_distr[1]:
              if gen not in strong_gens:
                  strong_gens.append(gen)
          return base, strong_gens
  
      def schreier_vector(self, alpha):
          """Computes the schreier vector for ``alpha``.
  
          The Schreier vector efficiently stores information
          about the orbit of ``alpha``. It can later be used to quickly obtain
          elements of the group that send ``alpha`` to a particular element
          in the orbit. Notice that the Schreier vector depends on the order
          in which the group generators are listed. For a definition, see [3].
          Since list indices start from zero, we adopt the convention to use
          "None" instead of 0 to signify that an element doesn't belong
          to the orbit.
          For the algorithm and its correctness, see [2], pp.78-80.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> from sympy.combinatorics.permutations import Permutation
          >>> a = Permutation([2, 4, 6, 3, 1, 5, 0])
          >>> b = Permutation([0, 1, 3, 5, 4, 6, 2])
          >>> G = PermutationGroup([a, b])
          >>> G.schreier_vector(0)
          [-1, None, 0, 1, None, 1, 0]
  
          See Also
          ========
  
          orbit
  
          """
          n = self.degree
          v = [None]*n
          v[alpha] = -1
          orb = [alpha]
          used = [False]*n
          used[alpha] = True
          gens = self.generators
          r = len(gens)
          for b in orb:
              for i in range(r):
                  temp = gens[i]._array_form[b]
                  if used[temp] is False:
                      orb.append(temp)
                      used[temp] = True
                      v[temp] = i
          return v
  
      def stabilizer(self, alpha):
          r"""Return the stabilizer subgroup of ``alpha``.
  
          The stabilizer of `\alpha` is the group `G_\alpha =
          \{g \in G | g(\alpha) = \alpha\}`.
          For a proof of correctness, see [1], p.79.
  
          Examples
          ========
  
          >>> from sympy.combinatorics import Permutation
          >>> Permutation.print_cyclic = True
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> from sympy.combinatorics.named_groups import DihedralGroup
          >>> G = DihedralGroup(6)
          >>> G.stabilizer(5)
          PermutationGroup([
              (5)(0 4)(1 3),
              (5)])
  
          See Also
          ========
  
          orbit
  
          """
          return PermGroup(_stabilizer(self._degree, self._generators, alpha))
  
      @property
      def strong_gens(self):
          """Return a strong generating set from the Schreier-Sims algorithm.
  
          A generating set `S = \{g_1, g_2, ..., g_t\}` for a permutation group
          `G` is a strong generating set relative to the sequence of points
          (referred to as a "base") `(b_1, b_2, ..., b_k)` if, for
          `1 \leq i \leq k` we have that the intersection of the pointwise
          stabilizer `G^{(i+1)} := G_{b_1, b_2, ..., b_i}` with `S` generates
          the pointwise stabilizer `G^{(i+1)}`. The concepts of a base and
          strong generating set and their applications are discussed in depth
          in [1], pp. 87-89 and [2], pp. 55-57.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.named_groups import DihedralGroup
          >>> D = DihedralGroup(4)
          >>> D.strong_gens
          [(0 1 2 3), (0 3)(1 2), (1 3)]
          >>> D.base
          [0, 1]
  
          See Also
          ========
  
          base, basic_transversals, basic_orbits, basic_stabilizers
  
          """
          if self._strong_gens == []:
              self.schreier_sims()
          return self._strong_gens
  
      def subgroup(self, gens):
          """
             Return the subgroup generated by `gens` which is a list of
             elements of the group
          """
  
          if not all([g in self for g in gens]):
              raise ValueError("The group doesn't contain the supplied generators")
  
          G = PermutationGroup(gens)
          return G
  
      def subgroup_search(self, prop, base=None, strong_gens=None, tests=None,
                          init_subgroup=None):
          """Find the subgroup of all elements satisfying the property ``prop``.
  
          This is done by a depth-first search with respect to base images that
          uses several tests to prune the search tree.
  
          Parameters
          ==========
  
          prop
              The property to be used. Has to be callable on group elements
              and always return ``True`` or ``False``. It is assumed that
              all group elements satisfying ``prop`` indeed form a subgroup.
          base
              A base for the supergroup.
          strong_gens
              A strong generating set for the supergroup.
          tests
              A list of callables of length equal to the length of ``base``.
              These are used to rule out group elements by partial base images,
              so that ``tests[l](g)`` returns False if the element ``g`` is known
              not to satisfy prop base on where g sends the first ``l + 1`` base
              points.
          init_subgroup
              if a subgroup of the sought group is
              known in advance, it can be passed to the function as this
              parameter.
  
          Returns
          =======
  
          res
              The subgroup of all elements satisfying ``prop``. The generating
              set for this group is guaranteed to be a strong generating set
              relative to the base ``base``.
  
          Examples
          ========
  
          >>> from sympy.combinatorics.named_groups import (SymmetricGroup,
          ... AlternatingGroup)
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> from sympy.combinatorics.testutil import _verify_bsgs
          >>> S = SymmetricGroup(7)
          >>> prop_even = lambda x: x.is_even
          >>> base, strong_gens = S.schreier_sims_incremental()
          >>> G = S.subgroup_search(prop_even, base=base, strong_gens=strong_gens)
          >>> G.is_subgroup(AlternatingGroup(7))
          True
          >>> _verify_bsgs(G, base, G.generators)
          True
  
          Notes
          =====
  
          This function is extremely lenghty and complicated and will require
          some careful attention. The implementation is described in
          [1], pp. 114-117, and the comments for the code here follow the lines
          of the pseudocode in the book for clarity.
  
          The complexity is exponential in general, since the search process by
          itself visits all members of the supergroup. However, there are a lot
          of tests which are used to prune the search tree, and users can define
          their own tests via the ``tests`` parameter, so in practice, and for
          some computations, it's not terrible.
  
          A crucial part in the procedure is the frequent base change performed
          (this is line 11 in the pseudocode) in order to obtain a new basic
          stabilizer. The book mentiones that this can be done by using
          ``.baseswap(...)``, however the current imlementation uses a more
          straightforward way to find the next basic stabilizer - calling the
          function ``.stabilizer(...)`` on the previous basic stabilizer.
  
          """
          # initialize BSGS and basic group properties
          def get_reps(orbits):
              # get the minimal element in the base ordering
              return [min(orbit, key = lambda x: base_ordering[x]) \
                for orbit in orbits]
  
          def update_nu(l):
              temp_index = len(basic_orbits[l]) + 1 -\
                           len(res_basic_orbits_init_base[l])
              # this corresponds to the element larger than all points
              if temp_index >= len(sorted_orbits[l]):
                  nu[l] = base_ordering[degree]
              else:
                  nu[l] = sorted_orbits[l][temp_index]
  
          if base is None:
              base, strong_gens = self.schreier_sims_incremental()
          base_len = len(base)
          degree = self.degree
          identity = _af_new(list(range(degree)))
          base_ordering = _base_ordering(base, degree)
          # add an element larger than all points
          base_ordering.append(degree)
          # add an element smaller than all points
          base_ordering.append(-1)
          # 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)
          # handle subgroup initialization and tests
          if init_subgroup is None:
              init_subgroup = PermutationGroup([identity])
          if tests is None:
              trivial_test = lambda x: True
              tests = []
              for i in range(base_len):
                  tests.append(trivial_test)
          # line 1: more initializations.
          res = init_subgroup
          f = base_len - 1
          l = base_len - 1
          # line 2: set the base for K to the base for G
          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_generators = res.generators
          res_basic_orbits_init_base = \
          [_orbit(degree, res_strong_gens_distr[i], res_base[i])\
           for i in range(base_len)]
          # initialize orbit representatives
          orbit_reps = [None]*base_len
          # line 4: orbit representatives for f-th basic stabilizer of K
          orbits = _orbits(degree, res_strong_gens_distr[f])
          orbit_reps[f] = get_reps(orbits)
          # line 5: remove the base point from the representatives to avoid
          # getting the identity element as a generator for K
          orbit_reps[f].remove(base[f])
          # line 6: more initializations
          c = [0]*base_len
          u = [identity]*base_len
          sorted_orbits = [None]*base_len
          for i in range(base_len):
              sorted_orbits[i] = basic_orbits[i][:]
              sorted_orbits[i].sort(key=lambda point: base_ordering[point])
          # line 7: initializations
          mu = [None]*base_len
          nu = [None]*base_len
          # this corresponds to the element smaller than all points
          mu[l] = degree + 1
          update_nu(l)
          # initialize computed words
          computed_words = [identity]*base_len
          # line 8: main loop
          while True:
              # apply all the tests
              while l < base_len - 1 and \
                  computed_words[l](base[l]) in orbit_reps[l] and \
                  base_ordering[mu[l]] < \
                  base_ordering[computed_words[l](base[l])] < \
                  base_ordering[nu[l]] and \
                      tests[l](computed_words):
                  # line 11: change the (partial) base of K
                  new_point = computed_words[l](base[l])
                  res_base[l] = new_point
                  new_stab_gens = _stabilizer(degree, res_strong_gens_distr[l],
                          new_point)
                  res_strong_gens_distr[l + 1] = new_stab_gens
                  # line 12: calculate minimal orbit representatives for the
                  # l+1-th basic stabilizer
                  orbits = _orbits(degree, new_stab_gens)
                  orbit_reps[l + 1] = get_reps(orbits)
                  # line 13: amend sorted orbits
                  l += 1
                  temp_orbit = [computed_words[l - 1](point) for point
                               in basic_orbits[l]]
                  temp_orbit.sort(key=lambda point: base_ordering[point])
                  sorted_orbits[l] = temp_orbit
                  # lines 14 and 15: update variables used minimality tests
                  new_mu = degree + 1
                  for i in range(l):
                      if base[l] in res_basic_orbits_init_base[i]:
                          candidate = computed_words[i](base[i])
                          if base_ordering[candidate] > base_ordering[new_mu]:
                              new_mu = candidate
                  mu[l] = new_mu
                  update_nu(l)
                  # line 16: determine the new transversal element
                  c[l] = 0
                  temp_point = sorted_orbits[l][c[l]]
                  gamma = computed_words[l - 1]._array_form.index(temp_point)
                  u[l] = transversals[l][gamma]
                  # update computed words
                  computed_words[l] = rmul(computed_words[l - 1], u[l])
              # lines 17 & 18: apply the tests to the group element found
              g = computed_words[l]
              temp_point = g(base[l])
              if l == base_len - 1 and \
                  base_ordering[mu[l]] < \
                  base_ordering[temp_point] < base_ordering[nu[l]] and \
                  temp_point in orbit_reps[l] and \
                  tests[l](computed_words) and \
                      prop(g):
                  # line 19: reset the base of K
                  res_generators.append(g)
                  res_base = base[:]
                  # line 20: recalculate basic orbits (and transversals)
                  res_strong_gens.append(g)
                  res_strong_gens_distr = _distribute_gens_by_base(res_base,
                                                            res_strong_gens)
                  res_basic_orbits_init_base = \
                  [_orbit(degree, res_strong_gens_distr[i], res_base[i]) \
                   for i in range(base_len)]
                  # line 21: recalculate orbit representatives
                  # line 22: reset the search depth
                  orbit_reps[f] = get_reps(orbits)
                  l = f
              # line 23: go up the tree until in the first branch not fully
              # searched
              while l >= 0 and c[l] == len(basic_orbits[l]) - 1:
                  l = l - 1
              # line 24: if the entire tree is traversed, return K
              if l == -1:
                  return PermutationGroup(res_generators)
              # lines 25-27: update orbit representatives
              if l < f:
                  # line 26
                  f = l
                  c[l] = 0
                  # line 27
                  temp_orbits = _orbits(degree, res_strong_gens_distr[f])
                  orbit_reps[f] = get_reps(temp_orbits)
                  # line 28: update variables used for minimality testing
                  mu[l] = degree + 1
                  temp_index = len(basic_orbits[l]) + 1 - \
                      len(res_basic_orbits_init_base[l])
                  if temp_index >= len(sorted_orbits[l]):
                      nu[l] = base_ordering[degree]
                  else:
                      nu[l] = sorted_orbits[l][temp_index]
              # line 29: set the next element from the current branch and update
              # accorndingly
              c[l] += 1
              if l == 0:
                  gamma  = sorted_orbits[l][c[l]]
              else:
                  gamma = computed_words[l - 1]._array_form.index(sorted_orbits[l][c[l]])
  
              u[l] = transversals[l][gamma]
              if l == 0:
                  computed_words[l] = u[l]
              else:
                  computed_words[l] = rmul(computed_words[l - 1], u[l])
  
      @property
      def transitivity_degree(self):
          """Compute the degree of transitivity of the group.
  
          A permutation group `G` acting on `\Omega = \{0, 1, ..., n-1\}` is
          ``k``-fold transitive, if, for any k points
          `(a_1, a_2, ..., a_k)\in\Omega` and any k points
          `(b_1, b_2, ..., b_k)\in\Omega` there exists `g\in G` such that
          `g(a_1)=b_1, g(a_2)=b_2, ..., g(a_k)=b_k`
          The degree of transitivity of `G` is the maximum ``k`` such that
          `G` is ``k``-fold transitive. ([8])
  
          Examples
          ========
  
          >>> from sympy.combinatorics.perm_groups import PermutationGroup
          >>> from sympy.combinatorics.permutations import Permutation
          >>> a = Permutation([1, 2, 0])
          >>> b = Permutation([1, 0, 2])
          >>> G = PermutationGroup([a, b])
          >>> G.transitivity_degree
          3
  
          See Also
          ========
          is_transitive, orbit
  
          """
          if self._transitivity_degree is None:
              n = self.degree
              G = self
              # if G is k-transitive, a tuple (a_0,..,a_k)
              # can be brought to (b_0,...,b_(k-1), b_k)
              # where b_0,...,b_(k-1) are fixed points;
              # consider the group G_k which stabilizes b_0,...,b_(k-1)
              # if G_k is transitive on the subset excluding b_0,...,b_(k-1)
              # then G is (k+1)-transitive
              for i in range(n):
                  orb = G.orbit((i))
                  if len(orb) != n - i:
                      self._transitivity_degree = i
                      return i
                  G = G.stabilizer(i)
              self._transitivity_degree = n
              return n
          else:
              return self._transitivity_degree
  
  
  def _orbit(degree, generators, alpha, action='tuples'):
      r"""Compute the orbit of alpha `\{g(\alpha) | g \in G\}` as a set.
  
      The time complexity of the algorithm used here is `O(|Orb|*r)` where
      `|Orb|` is the size of the orbit and ``r`` is the number of generators of
      the group. For a more detailed analysis, see [1], p.78, [2], pp. 19-21.
      Here alpha can be a single point, or a list of points.
  
      If alpha is a single point, the ordinary orbit is computed.
      if alpha is a list of points, there are three available options:
  
      'union' - computes the union of the orbits of the points in the list
      'tuples' - computes the orbit of the list interpreted as an ordered
      tuple under the group action ( i.e., g((1, 2, 3)) = (g(1), g(2), g(3)) )
      'sets' - computes the orbit of the list interpreted as a sets
  
      Examples
      ========
  
      >>> from sympy.combinatorics import Permutation
      >>> from sympy.combinatorics.perm_groups import PermutationGroup, _orbit
      >>> a = Permutation([1, 2, 0, 4, 5, 6, 3])
      >>> G = PermutationGroup([a])
      >>> _orbit(G.degree, G.generators, 0)
      {0, 1, 2}
      >>> _orbit(G.degree, G.generators, [0, 4], 'union')
      {0, 1, 2, 3, 4, 5, 6}
  
      See Also
      ========
  
      orbit, orbit_transversal
  
      """
      if not hasattr(alpha, '__getitem__'):
          alpha = [alpha]
  
      gens = [x._array_form for x in generators]
      if len(alpha) == 1 or action == 'union':
          orb = alpha
          used = [False]*degree
          for el in alpha:
              used[el] = True
          for b in orb:
              for gen in gens:
                  temp = gen[b]
                  if used[temp] == False:
                      orb.append(temp)
                      used[temp] = True
          return set(orb)
      elif action == 'tuples':
          alpha = tuple(alpha)
          orb = [alpha]
          used = {alpha}
          for b in orb:
              for gen in gens:
                  temp = tuple([gen[x] for x in b])
                  if temp not in used:
                      orb.append(temp)
                      used.add(temp)
          return set(orb)
      elif action == 'sets':
          alpha = frozenset(alpha)
          orb = [alpha]
          used = {alpha}
          for b in orb:
              for gen in gens:
                  temp = frozenset([gen[x] for x in b])
                  if temp not in used:
                      orb.append(temp)
                      used.add(temp)
          return {tuple(x) for x in orb}
  
  def _orbits(degree, generators):
      """Compute the orbits of G.
  
      If ``rep=False`` it returns a list of sets else it returns a list of
      representatives of the orbits
  
      Examples
      ========
  
      >>> from sympy.combinatorics.permutations import Permutation
      >>> from sympy.combinatorics.perm_groups import PermutationGroup, _orbits
      >>> a = Permutation([0, 2, 1])
      >>> b = Permutation([1, 0, 2])
      >>> _orbits(a.size, [a, b])
      [{0, 1, 2}]
      """
  
      seen = set()  # elements that have already appeared in orbits
      orbs = []
      sorted_I = list(range(degree))
      I = set(sorted_I)
      while I:
          i = sorted_I[0]
          orb = _orbit(degree, generators, i)
          orbs.append(orb)
          # remove all indices that are in this orbit
          I -= orb
          sorted_I = [i for i in sorted_I if i not in orb]
      return orbs
  
  def _orbit_transversal(degree, generators, alpha, pairs, af=False, slp=False):
      r"""Computes a transversal for the orbit of ``alpha`` as a set.
  
      generators   generators of the group ``G``
  
      For a permutation group ``G``, a transversal for the orbit
      `Orb = \{g(\alpha) | g \in G\}` is a set
      `\{g_\beta | g_\beta(\alpha) = \beta\}` for `\beta \in Orb`.
      Note that there may be more than one possible transversal.
      If ``pairs`` is set to ``True``, it returns the list of pairs
      `(\beta, g_\beta)`. For a proof of correctness, see [1], p.79
  
      if ``af`` is ``True``, the transversal elements are given in
      array form.
  
      If `slp` is `True`, a dictionary `{beta: slp_beta}` is returned
      for `\beta \in Orb` where `slp_beta` is a list of indices of the
      generators in `generators` s.t. if `slp_beta = [i_1 ... i_n]`
      `g_\beta = generators[i_n]*...*generators[i_1]`.
  
      Examples
      ========
  
      >>> from sympy.combinatorics import Permutation
      >>> Permutation.print_cyclic = True
      >>> from sympy.combinatorics.named_groups import DihedralGroup
      >>> from sympy.combinatorics.perm_groups import _orbit_transversal
      >>> G = DihedralGroup(6)
      >>> _orbit_transversal(G.degree, G.generators, 0, False)
      [(5), (0 1 2 3 4 5), (0 5)(1 4)(2 3), (0 2 4)(1 3 5), (5)(0 4)(1 3), (0 3)(1 4)(2 5)]
      """
  
      tr = [(alpha, list(range(degree)))]
      slp_dict = {alpha: []}
      used = [False]*degree
      used[alpha] = True
      gens = [x._array_form for x in generators]
      for x, px in tr:
          px_slp = slp_dict[x]
          for gen in gens:
              temp = gen[x]
              if used[temp] == False:
                  slp_dict[temp] = [gens.index(gen)] + px_slp
                  tr.append((temp, _af_rmul(gen, px)))
                  used[temp] = True
      if pairs:
          if not af:
              tr = [(x, _af_new(y)) for x, y in tr]
          if not slp:
              return tr
          return tr, slp_dict
  
      if af:
          tr = [y for _, y in tr]
          if not slp:
              return tr
          return tr, slp_dict
  
      tr = [_af_new(y) for _, y in tr]
      if not slp:
          return tr
      return tr, slp_dict
  
  def _stabilizer(degree, generators, alpha):
      r"""Return the stabilizer subgroup of ``alpha``.
  
      The stabilizer of `\alpha` is the group `G_\alpha =
      \{g \in G | g(\alpha) = \alpha\}`.
      For a proof of correctness, see [1], p.79.
  
      degree :       degree of G
      generators :   generators of G
  
      Examples
      ========
  
      >>> from sympy.combinatorics import Permutation
      >>> Permutation.print_cyclic = True
      >>> from sympy.combinatorics.perm_groups import _stabilizer
      >>> from sympy.combinatorics.named_groups import DihedralGroup
      >>> G = DihedralGroup(6)
      >>> _stabilizer(G.degree, G.generators, 5)
      [(5)(0 4)(1 3), (5)]
  
      See Also
      ========
  
      orbit
  
      """
      orb = [alpha]
      table = {alpha: list(range(degree))}
      table_inv = {alpha: list(range(degree))}
      used = [False]*degree
      used[alpha] = True
      gens = [x._array_form for x in generators]
      stab_gens = []
      for b in orb:
          for gen in gens:
              temp = gen[b]
              if used[temp] is False:
                  gen_temp = _af_rmul(gen, table[b])
                  orb.append(temp)
                  table[temp] = gen_temp
                  table_inv[temp] = _af_invert(gen_temp)
                  used[temp] = True
              else:
                  schreier_gen = _af_rmuln(table_inv[temp], gen, table[b])
                  if schreier_gen not in stab_gens:
                      stab_gens.append(schreier_gen)
      return [_af_new(x) for x in stab_gens]
  
  PermGroup = PermutationGroup
  