1from sympy.combinatorics.permutations import Permutation, _af_invert, _af_rmul
2from sympy.ntheory import isprime
3
4rmul = Permutation.rmul
5_af_new = Permutation._af_new
6
7############################################
8#
9# Utilities for computational group theory
10#
11############################################
12
13
14def _base_ordering(base, degree):
15    r"""
16    Order `\{0, 1, ..., n-1\}` so that base points come first and in order.
17
18    Parameters
19    ==========
20
21    ``base`` : the base
22    ``degree`` : the degree of the associated permutation group
23
24    Returns
25    =======
26
27    A list ``base_ordering`` such that ``base_ordering[point]`` is the
28    number of ``point`` in the ordering.
29
30    Examples
31    ========
32
33    >>> from sympy.combinatorics.named_groups import SymmetricGroup
34    >>> from sympy.combinatorics.util import _base_ordering
35    >>> S = SymmetricGroup(4)
36    >>> S.schreier_sims()
37    >>> _base_ordering(S.base, S.degree)
38    [0, 1, 2, 3]
39
40    Notes
41    =====
42
43    This is used in backtrack searches, when we define a relation `<<` on
44    the underlying set for a permutation group of degree `n`,
45    `\{0, 1, ..., n-1\}`, so that if `(b_1, b_2, ..., b_k)` is a base we
46    have `b_i << b_j` whenever `i<j` and `b_i << a` for all
47    `i\in\{1,2, ..., k\}` and `a` is not in the base. The idea is developed
48    and applied to backtracking algorithms in [1], pp.108-132. The points
49    that are not in the base are taken in increasing order.
50
51    References
52    ==========
53
54    .. [1] Holt, D., Eick, B., O'Brien, E.
55           "Handbook of computational group theory"
56
57    """
58    base_len = len(base)
59    ordering = [0]*degree
60    for i in range(base_len):
61        ordering[base[i]] = i
62    current = base_len
63    for i in range(degree):
64        if i not in base:
65            ordering[i] = current
66            current += 1
67    return ordering
68
69
70def _check_cycles_alt_sym(perm):
71    """
72    Checks for cycles of prime length p with n/2 < p < n-2.
73
74    Explanation
75    ===========
76
77    Here `n` is the degree of the permutation. This is a helper function for
78    the function is_alt_sym from sympy.combinatorics.perm_groups.
79
80    Examples
81    ========
82
83    >>> from sympy.combinatorics.util import _check_cycles_alt_sym
84    >>> from sympy.combinatorics.permutations import Permutation
85    >>> a = Permutation([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], [11, 12]])
86    >>> _check_cycles_alt_sym(a)
87    False
88    >>> b = Permutation([[0, 1, 2, 3, 4, 5, 6], [7, 8, 9, 10]])
89    >>> _check_cycles_alt_sym(b)
90    True
91
92    See Also
93    ========
94
95    sympy.combinatorics.perm_groups.PermutationGroup.is_alt_sym
96
97    """
98    n = perm.size
99    af = perm.array_form
100    current_len = 0
101    total_len = 0
102    used = set()
103    for i in range(n//2):
104        if not i in used and i < n//2 - total_len:
105            current_len = 1
106            used.add(i)
107            j = i
108            while af[j] != i:
109                current_len += 1
110                j = af[j]
111                used.add(j)
112            total_len += current_len
113            if current_len > n//2 and current_len < n - 2 and isprime(current_len):
114                return True
115    return False
116
117
118def _distribute_gens_by_base(base, gens):
119    r"""
120    Distribute the group elements ``gens`` by membership in basic stabilizers.
121
122    Explanation
123    ===========
124
125    Notice that for a base `(b_1, b_2, ..., b_k)`, the basic stabilizers
126    are defined as `G^{(i)} = G_{b_1, ..., b_{i-1}}` for
127    `i \in\{1, 2, ..., k\}`.
128
129    Parameters
130    ==========
131
132    ``base`` : a sequence of points in `\{0, 1, ..., n-1\}`
133    ``gens`` : a list of elements of a permutation group of degree `n`.
134
135    Returns
136    =======
137
138    List of length `k`, where `k` is
139    the length of ``base``. The `i`-th entry contains those elements in
140    ``gens`` which fix the first `i` elements of ``base`` (so that the
141    `0`-th entry is equal to ``gens`` itself). If no element fixes the first
142    `i` elements of ``base``, the `i`-th element is set to a list containing
143    the identity element.
144
145    Examples
146    ========
147
148    >>> from sympy.combinatorics.named_groups import DihedralGroup
149    >>> from sympy.combinatorics.util import _distribute_gens_by_base
150    >>> D = DihedralGroup(3)
151    >>> D.schreier_sims()
152    >>> D.strong_gens
153    [(0 1 2), (0 2), (1 2)]
154    >>> D.base
155    [0, 1]
156    >>> _distribute_gens_by_base(D.base, D.strong_gens)
157    [[(0 1 2), (0 2), (1 2)],
158     [(1 2)]]
159
160    See Also
161    ========
162
163    _strong_gens_from_distr, _orbits_transversals_from_bsgs,
164    _handle_precomputed_bsgs
165
166    """
167    base_len = len(base)
168    degree = gens[0].size
169    stabs = [[] for _ in range(base_len)]
170    max_stab_index = 0
171    for gen in gens:
172        j = 0
173        while j < base_len - 1 and gen._array_form[base[j]] == base[j]:
174            j += 1
175        if j > max_stab_index:
176            max_stab_index = j
177        for k in range(j + 1):
178            stabs[k].append(gen)
179    for i in range(max_stab_index + 1, base_len):
180        stabs[i].append(_af_new(list(range(degree))))
181    return stabs
182
183
184def _handle_precomputed_bsgs(base, strong_gens, transversals=None,
185                             basic_orbits=None, strong_gens_distr=None):
186    """
187    Calculate BSGS-related structures from those present.
188
189    Explanation
190    ===========
191
192    The base and strong generating set must be provided; if any of the
193    transversals, basic orbits or distributed strong generators are not
194    provided, they will be calculated from the base and strong generating set.
195
196    Parameters
197    ==========
198
199    ``base`` - the base
200    ``strong_gens`` - the strong generators
201    ``transversals`` - basic transversals
202    ``basic_orbits`` - basic orbits
203    ``strong_gens_distr`` - strong generators distributed by membership in basic
204    stabilizers
205
206    Returns
207    =======
208
209    ``(transversals, basic_orbits, strong_gens_distr)`` where ``transversals``
210    are the basic transversals, ``basic_orbits`` are the basic orbits, and
211    ``strong_gens_distr`` are the strong generators distributed by membership
212    in basic stabilizers.
213
214    Examples
215    ========
216
217    >>> from sympy.combinatorics.named_groups import DihedralGroup
218    >>> from sympy.combinatorics.util import _handle_precomputed_bsgs
219    >>> D = DihedralGroup(3)
220    >>> D.schreier_sims()
221    >>> _handle_precomputed_bsgs(D.base, D.strong_gens,
222    ... basic_orbits=D.basic_orbits)
223    ([{0: (2), 1: (0 1 2), 2: (0 2)}, {1: (2), 2: (1 2)}], [[0, 1, 2], [1, 2]], [[(0 1 2), (0 2), (1 2)], [(1 2)]])
224
225    See Also
226    ========
227
228    _orbits_transversals_from_bsgs, _distribute_gens_by_base
229
230    """
231    if strong_gens_distr is None:
232        strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
233    if transversals is None:
234        if basic_orbits is None:
235            basic_orbits, transversals = \
236                _orbits_transversals_from_bsgs(base, strong_gens_distr)
237        else:
238            transversals = \
239                _orbits_transversals_from_bsgs(base, strong_gens_distr,
240                                           transversals_only=True)
241    else:
242        if basic_orbits is None:
243            base_len = len(base)
244            basic_orbits = [None]*base_len
245            for i in range(base_len):
246                basic_orbits[i] = list(transversals[i].keys())
247    return transversals, basic_orbits, strong_gens_distr
248
249
250def _orbits_transversals_from_bsgs(base, strong_gens_distr,
251                                   transversals_only=False, slp=False):
252    """
253    Compute basic orbits and transversals from a base and strong generating set.
254
255    Explanation
256    ===========
257
258    The generators are provided as distributed across the basic stabilizers.
259    If the optional argument ``transversals_only`` is set to True, only the
260    transversals are returned.
261
262    Parameters
263    ==========
264
265    ``base`` - The base.
266    ``strong_gens_distr`` - Strong generators distributed by membership in basic
267    stabilizers.
268    ``transversals_only`` - bool
269        A flag switching between returning only the
270        transversals and both orbits and transversals.
271    ``slp`` -
272        If ``True``, return a list of dictionaries containing the
273        generator presentations of the elements of the transversals,
274        i.e. the list of indices of generators from ``strong_gens_distr[i]``
275        such that their product is the relevant transversal element.
276
277    Examples
278    ========
279
280    >>> from sympy.combinatorics.named_groups import SymmetricGroup
281    >>> from sympy.combinatorics.util import _distribute_gens_by_base
282    >>> S = SymmetricGroup(3)
283    >>> S.schreier_sims()
284    >>> strong_gens_distr = _distribute_gens_by_base(S.base, S.strong_gens)
285    >>> (S.base, strong_gens_distr)
286    ([0, 1], [[(0 1 2), (2)(0 1), (1 2)], [(1 2)]])
287
288    See Also
289    ========
290
291    _distribute_gens_by_base, _handle_precomputed_bsgs
292
293    """
294    from sympy.combinatorics.perm_groups import _orbit_transversal
295    base_len = len(base)
296    degree = strong_gens_distr[0][0].size
297    transversals = [None]*base_len
298    slps = [None]*base_len
299    if transversals_only is False:
300        basic_orbits = [None]*base_len
301    for i in range(base_len):
302        transversals[i], slps[i] = _orbit_transversal(degree, strong_gens_distr[i],
303                                 base[i], pairs=True, slp=True)
304        transversals[i] = dict(transversals[i])
305        if transversals_only is False:
306            basic_orbits[i] = list(transversals[i].keys())
307    if transversals_only:
308        return transversals
309    else:
310        if not slp:
311            return basic_orbits, transversals
312        return basic_orbits, transversals, slps
313
314
315def _remove_gens(base, strong_gens, basic_orbits=None, strong_gens_distr=None):
316    """
317    Remove redundant generators from a strong generating set.
318
319    Parameters
320    ==========
321
322    ``base`` - a base
323    ``strong_gens`` - a strong generating set relative to ``base``
324    ``basic_orbits`` - basic orbits
325    ``strong_gens_distr`` - strong generators distributed by membership in basic
326    stabilizers
327
328    Returns
329    =======
330
331    A strong generating set with respect to ``base`` which is a subset of
332    ``strong_gens``.
333
334    Examples
335    ========
336
337    >>> from sympy.combinatorics.named_groups import SymmetricGroup
338    >>> from sympy.combinatorics.util import _remove_gens
339    >>> from sympy.combinatorics.testutil import _verify_bsgs
340    >>> S = SymmetricGroup(15)
341    >>> base, strong_gens = S.schreier_sims_incremental()
342    >>> new_gens = _remove_gens(base, strong_gens)
343    >>> len(new_gens)
344    14
345    >>> _verify_bsgs(S, base, new_gens)
346    True
347
348    Notes
349    =====
350
351    This procedure is outlined in [1],p.95.
352
353    References
354    ==========
355
356    .. [1] Holt, D., Eick, B., O'Brien, E.
357           "Handbook of computational group theory"
358
359    """
360    from sympy.combinatorics.perm_groups import _orbit
361    base_len = len(base)
362    degree = strong_gens[0].size
363    if strong_gens_distr is None:
364        strong_gens_distr = _distribute_gens_by_base(base, strong_gens)
365    if basic_orbits is None:
366        basic_orbits = []
367        for i in range(base_len):
368            basic_orbit = _orbit(degree, strong_gens_distr[i], base[i])
369            basic_orbits.append(basic_orbit)
370    strong_gens_distr.append([])
371    res = strong_gens[:]
372    for i in range(base_len - 1, -1, -1):
373        gens_copy = strong_gens_distr[i][:]
374        for gen in strong_gens_distr[i]:
375            if gen not in strong_gens_distr[i + 1]:
376                temp_gens = gens_copy[:]
377                temp_gens.remove(gen)
378                if temp_gens == []:
379                    continue
380                temp_orbit = _orbit(degree, temp_gens, base[i])
381                if temp_orbit == basic_orbits[i]:
382                    gens_copy.remove(gen)
383                    res.remove(gen)
384    return res
385
386
387def _strip(g, base, orbits, transversals):
388    """
389    Attempt to decompose a permutation using a (possibly partial) BSGS
390    structure.
391
392    Explanation
393    ===========
394
395    This is done by treating the sequence ``base`` as an actual base, and
396    the orbits ``orbits`` and transversals ``transversals`` as basic orbits and
397    transversals relative to it.
398
399    This process is called "sifting". A sift is unsuccessful when a certain
400    orbit element is not found or when after the sift the decomposition
401    doesn't end with the identity element.
402
403    The argument ``transversals`` is a list of dictionaries that provides
404    transversal elements for the orbits ``orbits``.
405
406    Parameters
407    ==========
408
409    ``g`` - permutation to be decomposed
410    ``base`` - sequence of points
411    ``orbits`` - a list in which the ``i``-th entry is an orbit of ``base[i]``
412    under some subgroup of the pointwise stabilizer of `
413    `base[0], base[1], ..., base[i - 1]``. The groups themselves are implicit
414    in this function since the only information we need is encoded in the orbits
415    and transversals
416    ``transversals`` - a list of orbit transversals associated with the orbits
417    ``orbits``.
418
419    Examples
420    ========
421
422    >>> from sympy.combinatorics.named_groups import SymmetricGroup
423    >>> from sympy.combinatorics.permutations import Permutation
424    >>> from sympy.combinatorics.util import _strip
425    >>> S = SymmetricGroup(5)
426    >>> S.schreier_sims()
427    >>> g = Permutation([0, 2, 3, 1, 4])
428    >>> _strip(g, S.base, S.basic_orbits, S.basic_transversals)
429    ((4), 5)
430
431    Notes
432    =====
433
434    The algorithm is described in [1],pp.89-90. The reason for returning
435    both the current state of the element being decomposed and the level
436    at which the sifting ends is that they provide important information for
437    the randomized version of the Schreier-Sims algorithm.
438
439    References
440    ==========
441
442    .. [1] Holt, D., Eick, B., O'Brien, E."Handbook of computational group theory"
443
444    See Also
445    ========
446
447    sympy.combinatorics.perm_groups.PermutationGroup.schreier_sims
448    sympy.combinatorics.perm_groups.PermutationGroup.schreier_sims_random
449
450    """
451    h = g._array_form
452    base_len = len(base)
453    for i in range(base_len):
454        beta = h[base[i]]
455        if beta == base[i]:
456            continue
457        if beta not in orbits[i]:
458            return _af_new(h), i + 1
459        u = transversals[i][beta]._array_form
460        h = _af_rmul(_af_invert(u), h)
461    return _af_new(h), base_len + 1
462
463
464def _strip_af(h, base, orbits, transversals, j, slp=[], slps={}):
465    """
466    optimized _strip, with h, transversals and result in array form
467    if the stripped elements is the identity, it returns False, base_len + 1
468
469    j    h[base[i]] == base[i] for i <= j
470
471    """
472    base_len = len(base)
473    for i in range(j+1, base_len):
474        beta = h[base[i]]
475        if beta == base[i]:
476            continue
477        if beta not in orbits[i]:
478            if not slp:
479                return h, i + 1
480            return h, i + 1, slp
481        u = transversals[i][beta]
482        if h == u:
483            if not slp:
484                return False, base_len + 1
485            return False, base_len + 1, slp
486        h = _af_rmul(_af_invert(u), h)
487        if slp:
488            u_slp = slps[i][beta][:]
489            u_slp.reverse()
490            u_slp = [(i, (g,)) for g in u_slp]
491            slp = u_slp + slp
492    if not slp:
493        return h, base_len + 1
494    return h, base_len + 1, slp
495
496
497def _strong_gens_from_distr(strong_gens_distr):
498    """
499    Retrieve strong generating set from generators of basic stabilizers.
500
501    This is just the union of the generators of the first and second basic
502    stabilizers.
503
504    Parameters
505    ==========
506
507    ``strong_gens_distr`` - strong generators distributed by membership in basic
508    stabilizers
509
510    Examples
511    ========
512
513    >>> from sympy.combinatorics.named_groups import SymmetricGroup
514    >>> from sympy.combinatorics.util import (_strong_gens_from_distr,
515    ... _distribute_gens_by_base)
516    >>> S = SymmetricGroup(3)
517    >>> S.schreier_sims()
518    >>> S.strong_gens
519    [(0 1 2), (2)(0 1), (1 2)]
520    >>> strong_gens_distr = _distribute_gens_by_base(S.base, S.strong_gens)
521    >>> _strong_gens_from_distr(strong_gens_distr)
522    [(0 1 2), (2)(0 1), (1 2)]
523
524    See Also
525    ========
526
527    _distribute_gens_by_base
528
529    """
530    if len(strong_gens_distr) == 1:
531        return strong_gens_distr[0][:]
532    else:
533        result = strong_gens_distr[0]
534        for gen in strong_gens_distr[1]:
535            if gen not in result:
536                result.append(gen)
537        return result
538