# Source code for sympy.combinatorics.perm_groups

from random import randrange, choice
from math import log

from sympy.core import Basic
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

rmul = Permutation.rmul_with_af
_af_new = Permutation._af_new

[docs]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']

========

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

"""

def __new__(cls, *args, **kwargs):
"""The default constructor. Accepts Cycle and Permutation forms.
Removes duplicates unless dups keyword is False.
"""
args = list(args[0] if is_sequence(args[0]) else args)
if not args:
raise ValueError('must supply one or more permutations '
'to define the group')
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._basic_orbits = []
obj._transversals = []

# these attributes are assigned after running _random_pr_init
obj._random_gens = []
return obj

def __getitem__(self, i):
return self._generators[i]

def __len__(self):
return len(self._generators)

def __eq__(self, other):
"""Return True if self and other have the same generators.

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
return set(self.generators) == set(other.generators)

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([
Permutation(9)(0, 1, 2, 3, 4),
Permutation(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 = range(n1)
end = 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

========

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 = 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.

========

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.

========

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
[docs]    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]

========

strong_gens, basic_transversals, basic_orbits, basic_stabilizers

"""
if self._base == []:
self.schreier_sims()
return self._base

[docs]    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],
[Permutation(0, 1, 2, 3), Permutation(3)(0, 1), Permutation(1, 3, 2),
Permutation(2, 3), Permutation(1, 3)])

check that base, gens is a BSGS

>>> S1 = PermutationGroup(gens)
>>> _verify_bsgs(S1, base, gens)
True

========

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 = iter(Gamma).next()
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
[docs]    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]]

========

base, strong_gens, basic_transversals, basic_stabilizers

"""
if self._basic_orbits == []:
self.schreier_sims()
return self._basic_orbits

@property
[docs]    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([
Permutation(3)(0, 1, 2),
Permutation(1, 2, 3)])
PermutationGroup([
Permutation(1, 2, 3)])

========

base, strong_gens, basic_orbits, basic_transversals

"""

if self._transversals == []:
self.schreier_sims()
strong_gens = self._strong_gens
base = self._base
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
[docs]    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: Permutation(3),
1: Permutation(3)(0, 1, 2),
2: Permutation(3)(0, 2, 1),
3: Permutation(0, 3, 1)},
{1: Permutation(3),
2: Permutation(1, 2, 3),
3: Permutation(1, 3, 2)}]

========

strong_gens, base, basic_orbits, basic_stabilizers

"""

if self._transversals == []:
self.schreier_sims()
return self._transversals

[docs]    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

========

centralizer

Notes
=====

This is a naive implementation that is a straightforward application
of .centralizer()

"""
return self.centralizer(self)

[docs]    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

========

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(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]))

[docs]    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

========

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

[docs]    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 = 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

[docs]    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)
Permutation(7)(2, 4)(3, 5)

========

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

[docs]    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
[docs]    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())
[Permutation(2), Permutation(2)(0, 1)]

========

order
"""
return self._degree

[docs]    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

========

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

[docs]    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]]

========

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 = range(degree)
for i in range(r):
for j in range(r):
p1 = gens[i]
p2 = gens[j]
c = range(degree)
for k in rng:
c[p2[p1[k]]] = p1[p2[k]]
ct = tuple(c)
if not ct in set_commutators:
cms = [_af_new(p) for p in set_commutators]
G2 = self.normal_closure(cms)
return G2

[docs]    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 not
true groups:

>>> G = tetrahedron.pgroup
>>> G.is_group()
False

But 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([
Permutation(0, 1)(2, 3),
Permutation(3),
Permutation(1, 2, 3),
Permutation(1, 3, 2),
Permutation(0, 3, 1),
Permutation(0, 2, 3),
Permutation(0, 3)(1, 2),
Permutation(0, 1, 3),
Permutation(3)(0, 2, 1),
Permutation(0, 3, 2),
Permutation(3)(0, 1, 2),
Permutation(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)

[docs]    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 = range(self.degree)
order = 0
element_list = [idn]
set_element_list = set([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)
N.append(ap)
self._order = len(element_list)

[docs]    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
raise StopIteration
if len(u) == 1:
for i in basic_orbits[0]:
if af:
yield u[0][i]._array_form
else:
yield u[0][i]
raise StopIteration

u = list(reversed(u))
basic_orbits = basic_orbits[::-1]
# stg stack of group elements
stg = [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]:
#count_b += 1
if h == 0:
raise StopIteration
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
[docs]    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
[Permutation(1, 2), Permutation(2)(0, 1)]

"""
return self._generators

[docs]    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

========

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
[docs]    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

[docs]    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

========

_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
[docs]    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

========

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(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

[docs]    def is_normal(self, gr):
"""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.generator and
g2 g2 in G.generator

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

"""
gens2 = [p._array_form for p in 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 self.coset_factor(p, True):
return False
return True

[docs]    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

========

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
[docs]    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

========

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(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

[docs]    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)

[docs]    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
[docs]    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

[docs]    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

========

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
[docs]    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

========

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

[docs]    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]

========

_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 = 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:
temp = 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(temp, parents)
# union has side effects: performs union by rank on the list
# of representatives
temp = self._union_find_merge(gen(temp), 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

[docs]    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

========

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(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]))

[docs]    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)
set([0, 1, 2])
>>> G.orbit([0,4], 'union')
set([0, 1, 2, 3, 4, 5, 6])

========

orbit_transversal

"""
return _orbit(self.degree, self.generators, alpha, action)

[docs]    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)
Permutation(0, 4, 1, 2, 3)

========

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(range(self._degree))

[docs]    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)
[Permutation(5),
Permutation(0, 1, 2, 3, 4, 5),
Permutation(0, 5)(1, 4)(2, 3),
Permutation(0, 2, 4)(1, 3, 5),
Permutation(5)(0, 4)(1, 3),
Permutation(0, 3)(1, 4)(2, 5)]

========

orbit

"""
return _orbit_transversal(self._degree, self.generators, alpha, pairs)

[docs]    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()
[set([0, 2, 3, 4, 6]), set([1, 5])]
"""
return _orbits(self._degree, self._generators)

[docs]    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())
[Permutation(2), Permutation(2)(0, 1)]

>>> a = Permutation([0, 2, 1])
>>> b = Permutation([1, 0, 2])
>>> G = PermutationGroup([a, b])
>>> G.order()
6

========

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

[docs]    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

========

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(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)

[docs]    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])
Permutation(0, 1)(2, 3)
>>> G.make_perm(3, [0, 1, 0])
Permutation(0, 2, 3, 1)
>>> G.make_perm([0, 1, 0])
Permutation(0, 2, 3, 1)

========

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)

result = Permutation(range(self.degree))
m = len(self)
for i in range(n):
p = self[randrange(m)]
result = rmul(result, p)
return result

[docs]    def random(self, af=False):
"""Return a random group element
"""
rank = randrange(self.order())
return self.coset_unrank(rank, af)

[docs]    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.

========

_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])

[docs]    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

========

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)

[docs]    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: Permutation(2)(0, 1), 1: Permutation(2), 2: Permutation(1, 2)},
{0: Permutation(2), 2: Permutation(0, 2)}]
"""
if self._transversals:
return
base, strong_gens = self.schreier_sims_incremental()
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 = _orbits_transversals_from_bsgs(base,\
strong_gens_distr)
self._transversals = transversals
self._basic_orbits = [sorted(x) for x in basic_orbits]

[docs]    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.

========

schreier_sims, schreier_sims_random

"""
if base is None:
base = []
if gens is None:
gens = self.generators[:]
degree = self.degree
id_af = range(degree)
# handle the trivial group
if len(gens) == 1 and gens[0].is_Identity:
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)
# initialize the basic stabilizers, basic orbits and basic transversals
orbs = {}
transversals = {}
base_len = len(_base)
for i in range(base_len):
transversals[i] = dict(_orbit_transversal(degree, strong_gens_distr[i],
_base[i], pairs=True, af=True))
orbs[i] = 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 transversals[i].items():
for gen in strong_gens_distr[i]:
gb = gen._array_form[beta]
u1 = transversals[i][gb]
g1 = _af_rmul(gen._array_form, u_beta)
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)
h, j = _strip_af(schreier_gen, _base, orbs, transversals, i)
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)
for l in range(i + 1, j):
strong_gens_distr[l].append(h)
transversals[l] =\
dict(_orbit_transversal(degree, strong_gens_distr[l],
_base[l], pairs=True, af=True))
orbs[l] = 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 = []
for gens in strong_gens_distr:
for gen in gens:
if gen not in strong_gens:
strong_gens.append(gen)
return _base, strong_gens

[docs]    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}.

========

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] = 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] = 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

[docs]    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]

========

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

[docs]    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([
Permutation(5)(0, 4)(1, 3),
Permutation(5)])

========

orbit

"""
return PermGroup(_stabilizer(self._degree, self._generators, alpha))

@property
[docs]    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
[Permutation(0, 1, 2, 3), Permutation(0, 3)(1, 2), Permutation(1, 3)]
>>> D.base
[0, 1]

========

base, basic_transversals, basic_orbits, basic_stabilizers

"""
if self._strong_gens == []:
self.schreier_sims()
return self._strong_gens

[docs]    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(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
[docs]    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

========
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

[docs]    def is_group(self):
"""Return True if the group meets three criteria: identity is present,
the inverse of every element is also an element, and the product of
any two elements is also an element. If any of the tests fail, False
is returned.

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 not
a true group:

>>> G = tetrahedron.pgroup
>>> G.is_group()
False

But the group generated by the permutations in the tetrahedron
pgroup is a proper group:

>>> H = PermutationGroup(list(G.generate()))
>>> H.is_group()
True

The identity permutation is present:

>>> H.has(Permutation(G.degree - 1))
True

The product of any two elements from the group is also in the group:

>>> from sympy import TableForm
>>> g = list(H)
>>> n = len(g)
>>> m = []
>>> for i in g:
...     m.append([g.index(i*H) for H in g])
...
| 0  1  2  3  4  5  6  7  8  9  10 11
----------------------------------------
0 | 11 0  8  10 6  2  7  4  5  3  9  1
1 | 0  1  2  3  4  5  6  7  8  9  10 11
2 | 6  2  7  4  5  3  9  1  11 0  8  10
3 | 5  3  9  1  11 0  8  10 6  2  7  4
4 | 3  4  0  2  10 6  11 8  9  7  1  5
5 | 4  5  6  7  8  9  10 11 0  1  2  3
6 | 10 6  11 8  9  7  1  5  3  4  0  2
7 | 9  7  1  5  3  4  0  2  10 6  11 8
8 | 7  8  4  6  2  10 3  0  1  11 5  9
9 | 8  9  10 11 0  1  2  3  4  5  6  7
10 | 2  10 3  0  1  11 5  9  7  8  4  6
11 | 1  11 5  9  7  8  4  6  2  10 3  0
>>>
The entries in the table give the element in the group corresponding
to the product of a given column element and row element:

>>> g[3]*g[2] == g[9]
True

The inverse of every element is also in the group:

>>> TableForm([[g.index(~gi) for gi in g]], headings=[[], range(n)],
...    wipe_zeros=False)
0  1 2 3 4  5 6 7 8 9 10 11
---------------------------
11 1 7 3 10 9 6 2 8 5 4  0

So we see that g[1] and g[3] are equivalent to their inverse while
g[7] == ~g[2].
"""
# identity present
I = Permutation(size=self.degree)
for g in self:
if g == I:
break
else:
return False

# associativity already holds: a*(b*c) == (a*b)*c for permutations

# inverse of each is present
if not all(self.has(~a) for a in self):
return False

# closure
for a in self:
for b in self:
if not self.has(a*b):
return False
return True

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)
set([0, 1, 2])
>>> _orbit(G.degree, G.generators, [0, 4], 'union')
set([0, 1, 2, 3, 4, 5, 6])

========

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 = set([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)
return set(orb)
elif action == 'sets':
alpha = frozenset(alpha)
orb = [alpha]
used = set([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)
return set([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])
[set([0, 1, 2])]
"""

seen = set()  # elements that have already appeared in orbits
orbs = []
sorted_I = 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):
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

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)
[Permutation(5),
Permutation(0, 1, 2, 3, 4, 5),
Permutation(0, 5)(1, 4)(2, 3),
Permutation(0, 2, 4)(1, 3, 5),
Permutation(5)(0, 4)(1, 3),
Permutation(0, 3)(1, 4)(2, 5)]
"""

tr = [(alpha, range(degree))]
used = [False]*degree
used[alpha] = True
gens = [x._array_form for x in generators]
for x, px in tr:
for gen in gens:
temp = gen[x]
if used[temp] == False:
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]
return tr

if af:
return [y for _, y in tr]

return [_af_new(y) for _, y in tr]

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)
[Permutation(5)(0, 4)(1, 3), Permutation(5)]

========

orbit

"""
orb = [alpha]
table = {alpha: range(degree)}
table_inv = {alpha: 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