def _strip_af(h, base, orbits, transversals, j, slp=[], slps={}):
    """
    optimized _strip, with h, transversals and result in array form
    if the stripped elements is the identity, it returns False, base_len + 1

    j    h[base[i]] == base[i] for i <= j
    """
    base_len = len(base)
    for i in range(j+1, base_len):
        beta = h[base[i]]
        if beta == base[i]:
            continue
        if beta not in orbits[i]:
            if not slp:
                return h, i + 1
            return h, i + 1, slp
        u = transversals[i][beta]
        if h == u:
            if not slp:
                return False, base_len + 1
            return False, base_len + 1, slp
        h = _af_rmul(_af_invert(u), h)
        if slp:
            u_slp = slps[i][beta][:]
            u_slp.reverse()
            u_slp = [(i, (g,)) for g in u_slp]
            slp = u_slp + slp
    if not slp:
        return h, base_len + 1
    return h, base_len + 1, slp

def test_mul():
    a, b = [0, 2, 1, 3], [0, 1, 3, 2]
    assert _af_rmul(a, b) == [0, 2, 3, 1]
    assert _af_rmuln(a, b, range(4)) == [0, 2, 3, 1]
    assert rmul(Permutation(a), Permutation(b)).array_form == [0, 2, 3, 1]

    a = Permutation([0, 2, 1, 3])
    b = (0, 1, 3, 2)
    c = (3, 1, 2, 0)
    assert Permutation.rmul(a, b, c) == Permutation([1, 2, 3, 0])
    assert Permutation.rmul(a, c) == Permutation([3, 2, 1, 0])
    raises(TypeError, lambda: Permutation.rmul(b, c))

    n = 6
    m = 8
    a = [Permutation.unrank_nonlex(n, i).array_form for i in range(m)]
    h = range(n)
    for i in range(m):
        h = _af_rmul(h, a[i])
        h2 = _af_rmuln(*a[:i + 1])
        assert h == h2

def _strip_af(h, base, orbits, transversals, j):
    """
    optimized _strip, with h, transversals and result in array form
    if the stripped elements is the identity, it returns False, base_len + 1

    j    h[base[i]] == base[i] for i <= j
    """
    base_len = len(base)
    for i in range(j+1, base_len):
        beta = h[base[i]]
        if beta == base[i]:
            continue
        if beta not in orbits[i]:
            return h, i + 1
        u = transversals[i][beta]
        if h == u:
            return False, base_len + 1
        h = _af_rmul(_af_invert(u), h)
    return h, base_len + 1

def test_Permutation():
    # don't auto fill 0
    raises(ValueError, lambda: Permutation([1]))
    p = Permutation([0, 1, 2, 3])
    # call as bijective
    assert [p(i) for i in range(p.size)] == list(p)
    # call as operator
    assert p(range(p.size)) == list(p)
    # call as function
    assert list(p(1, 2)) == [0, 2, 1, 3]
    # conversion to list
    assert list(p) == range(4)
    # cycle form with size
    assert Permutation([[1, 2]], size=4) == Permutation([[1, 2], [0], [3]])
    # random generation
    assert Permutation.random(2) in (Permutation([1, 0]), Permutation([0, 1]))

    p = Permutation([2, 5, 1, 6, 3, 0, 4])
    q = Permutation([[1], [0, 3, 5, 6, 2, 4]])
    assert len(set([p, p])) == 1
    r = Permutation([1, 3, 2, 0, 4, 6, 5])
    ans = Permutation(_af_rmuln(*[w.array_form for w in (p, q, r)])).array_form
    assert rmul(p, q, r).array_form == ans
    # make sure no other permutation of p, q, r could have given
    # that answer
    for a, b, c in permutations((p, q, r)):
        if (a, b, c) == (p, q, r):
            continue
        assert rmul(a, b, c).array_form != ans

    assert p.support() == range(7)
    assert q.support() == [0, 2, 3, 4, 5, 6]
    assert Permutation(p.cyclic_form).array_form == p.array_form
    assert p.cardinality == 5040
    assert q.cardinality == 5040
    assert q.cycles == 2
    assert rmul(q, p) == Permutation([4, 6, 1, 2, 5, 3, 0])
    assert rmul(p, q) == Permutation([6, 5, 3, 0, 2, 4, 1])
    assert _af_rmul(p.array_form, q.array_form) == \
        [6, 5, 3, 0, 2, 4, 1]

    assert rmul(Permutation([[1, 2, 3], [0, 4]]),
                Permutation([[1, 2, 4], [0], [3]])).cyclic_form == \
        [[0, 4, 2], [1, 3]]
    assert q.array_form == [3, 1, 4, 5, 0, 6, 2]
    assert q.cyclic_form == [[0, 3, 5, 6, 2, 4]]
    assert q.full_cyclic_form == [[0, 3, 5, 6, 2, 4], [1]]
    assert p.cyclic_form == [[0, 2, 1, 5], [3, 6, 4]]
    t = p.transpositions()
    assert t == [(0, 5), (0, 1), (0, 2), (3, 4), (3, 6)]
    assert Permutation.rmul(*[Permutation(Cycle(*ti)) for ti in (t)])
    assert Permutation([1, 0]).transpositions() == [(0, 1)]

    assert p**13 == p
    assert q**0 == Permutation(range(q.size))
    assert q**-2 == ~q**2
    assert q**2 == Permutation([5, 1, 0, 6, 3, 2, 4])
    assert q**3 == q**2*q
    assert q**4 == q**2*q**2

    a = Permutation(1, 3)
    b = Permutation(2, 0, 3)
    I = Permutation(3)
    assert ~a == a**-1
    assert a*~a == I
    assert a*b**-1 == a*~b

    ans = Permutation(0, 5, 3, 1, 6)(2, 4)
    assert (p + q.rank()).rank() == ans.rank()
    assert (p + q.rank())._rank == ans.rank()
    assert (q + p.rank()).rank() == ans.rank()
    raises(TypeError, lambda: p + Permutation(range(10)))

    assert (p - q.rank()).rank() == Permutation(0, 6, 3, 1, 2, 5, 4).rank()
    assert p.rank() - q.rank() < 0  # for coverage: make sure mod is used
    assert (q - p.rank()).rank() == Permutation(1, 4, 6, 2)(3, 5).rank()

    assert p*q == Permutation(_af_rmuln(*[list(w) for w in (q, p)]))
    assert p*Permutation([]) == p
    assert Permutation([])*p == p
    assert p*Permutation([[0, 1]]) == Permutation([2, 5, 0, 6, 3, 1, 4])
    assert Permutation([[0, 1]])*p == Permutation([5, 2, 1, 6, 3, 0, 4])

    pq = p^q
    assert pq == Permutation([5, 6, 0, 4, 1, 2, 3])
    assert pq == rmul(q, p, ~q)
    qp = q^p
    assert qp == Permutation([4, 3, 6, 2, 1, 5, 0])
    assert qp == rmul(p, q, ~p)
    raises(ValueError, lambda: p^Permutation([]))

    assert p.commutator(q) == Permutation(0, 1, 3, 4, 6, 5, 2)
    assert q.commutator(p) == Permutation(0, 2, 5, 6, 4, 3, 1)
    assert p.commutator(q) == ~q.commutator(p)
    raises(ValueError, lambda: p.commutator(Permutation([])))

    assert len(p.atoms()) == 7
    assert q.atoms() == set([0, 1, 2, 3, 4, 5, 6])

    assert p.inversion_vector() == [2, 4, 1, 3, 1, 0]
#.........这里部分代码省略.........

def _strip(g, base, orbits, transversals):
    """
    Attempt to decompose a permutation using a (possibly partial) BSGS
    structure.

    This is done by treating the sequence ``base`` as an actual base, and
    the orbits ``orbits`` and transversals ``transversals`` as basic orbits and
    transversals relative to it.

    This process is called "sifting". A sift is unsuccessful when a certain
    orbit element is not found or when after the sift the decomposition
    doesn't end with the identity element.

    The argument ``transversals`` is a list of dictionaries that provides
    transversal elements for the orbits ``orbits``.

    Parameters
    ==========

    ``g`` - permutation to be decomposed
    ``base`` - sequence of points
    ``orbits`` - a list in which the ``i``-th entry is an orbit of ``base[i]``
    under some subgroup of the pointwise stabilizer of `
    `base[0], base[1], ..., base[i - 1]``. The groups themselves are implicit
    in this function since the only infromation we need is encoded in the orbits
    and transversals
    ``transversals`` - a list of orbit transversals associated with the orbits
    ``orbits``.

    Examples
    ========

    >>> from sympy.combinatorics import Permutation
    >>> Permutation.print_cyclic = True
    >>> from sympy.combinatorics.named_groups import SymmetricGroup
    >>> from sympy.combinatorics.permutations import Permutation
    >>> from sympy.combinatorics.util import _strip
    >>> S = SymmetricGroup(5)
    >>> S.schreier_sims()
    >>> g = Permutation([0, 2, 3, 1, 4])
    >>> _strip(g, S.base, S.basic_orbits, S.basic_transversals)
    (Permutation(4), 5)

    Notes
    =====

    The algorithm is described in [1],pp.89-90. The reason for returning
    both the current state of the element being decomposed and the level
    at which the sifting ends is that they provide important information for
    the randomized version of the Schreier-Sims algorithm.

    References
    ==========

    [1] Holt, D., Eick, B., O'Brien, E.
    "Handbook of computational group theory"

    See Also
    ========

    sympy.combinatorics.perm_groups.PermutationGroup.schreier_sims

    sympy.combinatorics.perm_groups.PermutationGroup.schreier_sims_random

    """
    h = g.array_form
    base_len = len(base)
    for i in range(base_len):
        beta = h[base[i]]
        if beta == base[i]:
            continue
        if beta not in orbits[i]:
            return _af_new(h), i + 1
        u = transversals[i][beta].array_form
        h = _af_rmul(_af_invert(u), h)
    return _af_new(h), base_len + 1

def canonical_free(base, gens, g, num_free):
    """
    canonicalization of a tensor with respect to free indices
    choosing the minimum with respect to lexicographical ordering
    in the free indices

    ``base``, ``gens``  BSGS for slot permutation group
    ``g``               permutation representing the tensor
    ``num_free``        number of free indices
    The indices must be ordered with first the free indices

    see explanation in double_coset_can_rep
    The algorithm is a variation of the one given in [2].

    Examples
    ========

    >>> from sympy.combinatorics import Permutation
    >>> from sympy.combinatorics.tensor_can import canonical_free
    >>> gens = [[1,0,2,3,5,4], [2,3,0,1,4,5],[0,1,3,2,5,4]]
    >>> gens = [Permutation(h) for h in gens]
    >>> base = [0, 2]
    >>> g = Permutation([2, 1, 0, 3, 4, 5])
    >>> canonical_free(base, gens, g, 4)
    [0, 3, 1, 2, 5, 4]

    Consider the product of Riemann tensors
    ``T = R^{a}_{d0}^{d1,d2}*R_{d2,d1}^{d0,b}``
    The order of the indices is ``[a,b,d0,-d0,d1,-d1,d2,-d2]``
    The permutation corresponding to the tensor is
    ``g = [0,3,4,6,7,5,2,1,8,9]``

    In particular ``a`` is position ``0``, ``b`` is in position ``9``.
    Use the slot symmetries to get `T` is a form which is the minimal
    in lexicographic order in the free indices ``a`` and ``b``, e.g.
    ``-R^{a}_{d0}^{d1,d2}*R^{b,d0}_{d2,d1}`` corresponding to
    ``[0, 3, 4, 6, 1, 2, 7, 5, 9, 8]``

    >>> from sympy.combinatorics.tensor_can import riemann_bsgs, tensor_gens
    >>> base, gens = riemann_bsgs
    >>> size, sbase, sgens = tensor_gens(base, gens, [[],[]], 0)
    >>> g = Permutation([0,3,4,6,7,5,2,1,8,9])
    >>> canonical_free(sbase, [Permutation(h) for h in sgens], g, 2)
    [0, 3, 4, 6, 1, 2, 7, 5, 9, 8]
    """
    g = g.array_form
    size = len(g)
    if not base:
        return g[:]

    transversals = get_transversals(base, gens)
    m = len(base)
    for x in sorted(g[:-2]):
        if x not in base:
            base.append(x)
    h = g
    for i, transv in enumerate(transversals):
        b = base[i]
        h_i = [size]*num_free
        # find the element s in transversals[i] such that
        # _af_rmul(h, s) has its free elements with the lowest position in h
        s = None
        for sk in transv.values():
            h1 = _af_rmul(h, sk)
            hi = [h1.index(ix) for ix in range(num_free)]
            if hi < h_i:
                h_i = hi
                s = sk
        if s:
            h = _af_rmul(h, s)
    return h

#.........这里部分代码省略.........

          `T^{a0 a1 a2 a3 a4 a5} = -T^{a4 a1 a2 a3 a0 a5}`

    and symmetric metric, find the tensor equivalent to it which
    is the lowest under the ordering of indices:
    lexicographic ordering `d1, d2, d3` then and contravariant index
    before covariant index; that is the canonical form of the tensor.

    The canonical form is `-T^{d1 d2 d3}{}_{d1 d2 d3}`
    obtained using `T^{a0 a1 a2 a3 a4 a5} = -T^{a2 a1 a0 a3 a4 a5}`.

    To convert this problem in the input for this function,
    use the following labelling of the index names
    (- for covariant for short) `d1, -d1, d2, -d2, d3, -d3`

    `T^{d3 d2 d1}{}_{d1 d2 d3}` corresponds to `g = [4,2,0,1,3,5,6,7]`
    where the last two indices are for the sign

    `sgens = [Permutation(0,2)(6,7), Permutation(0,4)(6,7)]`

    sgens[0] is the slot symmetry `-(0,2)`
    `T^{a0 a1 a2 a3 a4 a5} = -T^{a2 a1 a0 a3 a4 a5}`

    sgens[1] is the slot symmetry `-(0,4)`
    `T^{a0 a1 a2 a3 a4 a5} = -T^{a4 a1 a2 a3 a0 a5}`

    The dummy symmetry group D is generated by the strong base generators
    `[(0,1),(2,3),(4,5),(0,1)(2,3),(2,3)(4,5)]`

    The dummy symmetry acts from the left
    `d = [1,0,2,3,4,5,6,7]`  exchange `d1 -> -d1`
    `T^{d3 d2 d1}{}_{d1 d2 d3} == T^{d3 d2}{}_{d1}{}^{d1}{}_{d2 d3}`

    `g=[4,2,0,1,3,5,6,7]  -> [4,2,1,0,3,5,6,7] = _af_rmul(d, g)`
    which differs from `_af_rmul(g, d)`.

    The slot symmetry acts from the right
    `s = [2,1,0,3,4,5,7,6]`  exchanges slots 0 and 2 and changes sign
    `T^{d3 d2 d1}{}_{d1 d2 d3} == -T^{d1 d2 d3}{}_{d1 d2 d3}`

    `g=[4,2,0,1,3,5,6,7]  -> [0,2,4,1,3,5,7,6] = _af_rmul(g, s)`

    Example in which the tensor is zero, same slot symmetries as above:
    `T^{d3}{}_{d1,d2}{}^{d1}{}_{d3}{}^{d2}`

    `= -T^{d3}{}_{d1,d3}{}^{d1}{}_{d2}{}^{d2}`   under slot symmetry `-(2,4)`;

    `= T_{d3 d1}{}^{d3}{}^{d1}{}_{d2}{}^{d2}`    under slot symmetry `-(0,2)`;

    `= T^{d3}{}_{d1 d3}{}^{d1}{}_{d2}{}^{d2}`    symmetric metric;

    `= 0`  since two of these lines have tensors differ only for the sign.

    The double coset D*g*S consists of permutations `h = d*g*s` corresponding
    to equivalent tensors; if there are two `h` which are the same apart
    from the sign, return zero; otherwise
    choose as representative the tensor with indices
    ordered lexicographically according to `[d1, -d1, d2, -d2, d3, -d3]`
    that is `rep = min(D*g*S) = min([d*g*s for d in D for s in S])`

    The indices are fixed one by one; first choose the lowest index
    for slot 0, then the lowest remaining index for slot 1, etc.
    Doing this one obtains a chain of stabilizers

    `S -> S_{b0} -> S_{b0,b1} -> ...` and
    `D -> D_{p0} -> D_{p0,p1} -> ...`

def canonicalize_naive(g, dummies, sym, *v):
    """
    canonicalize tensor formed by tensors of the different types

    g  permutation representing the tensor
    dummies  list of dummy indices
    msym symmetry of the metric

    v is a list of (base_i, gens_i, n_i, sym_i) for tensors of type `i`
    base_i, gens_i BSGS for tensors of this type
    n_i  number ot tensors of type `i`

    sym_i symmetry under exchange of two component tensors of type `i`
          None  no symmetry
          0     commuting
          1     anticommuting

    Return 0 if the tensor is zero, else return the array form of
    the permutation representing the canonical form of the tensor.

    Examples
    ========

    >>> from sympy.combinatorics.testutil import canonicalize_naive
    >>> from sympy.combinatorics.tensor_can import get_symmetric_group_sgs
    >>> from sympy.combinatorics import Permutation, PermutationGroup
    >>> g = Permutation([1,3,2,0,4,5])
    >>> base2, gens2 = get_symmetric_group_sgs(2)
    >>> canonicalize_naive(g, [2, 3], 0, (base2, gens2, 2, 0))
    [0, 2, 1, 3, 4, 5]
    """
    from sympy.combinatorics.perm_groups import PermutationGroup
    from sympy.combinatorics.tensor_can import gens_products, dummy_sgs
    from sympy.combinatorics.permutations import Permutation, _af_rmul
    v1 = []
    for i in range(len(v)):
        base_i, gens_i, n_i, sym_i = v[i]
        v1.append((base_i, gens_i, [[]]*n_i, sym_i))
    size, sbase, sgens = gens_products(*v1)
    dgens = dummy_sgs(dummies, sym, size-2)
    if isinstance(sym, int):
        num_types = 1
        dummies = [dummies]
        sym = [sym]
    else:
        num_types = len(sym)
    dgens = []
    for i in range(num_types):
        dgens.extend(dummy_sgs(dummies[i], sym[i], size - 2))
    S = PermutationGroup(sgens)
    D = PermutationGroup([Permutation(x) for x in dgens])
    dlist = list(D.generate(af=True))
    g = g.array_form
    st = set()
    for s in S.generate(af=True):
        h = _af_rmul(g, s)
        for d in dlist:
            q = tuple(_af_rmul(d, h))
            st.add(q)
    a = list(st)
    a.sort()
    prev = (0,)*size
    for h in a:
        if h[:-2] == prev[:-2]:
            if h[-1] != prev[-1]:
                return 0
        prev = h
    return list(a[0])

def canonical_free(base, gens, g, num_free):
    """
    canonicalization of a tensor with respect to free indices
    choosing the minimum with respect to lexicographical ordering

    base, gens  BSGS for slot permutation group
    g           permutation representing the tensor
    num_free    number of free indices
    The indices must be ordered with first the free indices

    see explanation in double_coset_can_rep
    The algorithm is given in [2]

    Examples
    ========
    >>> from sympy.combinatorics import Permutation
    >>> from sympy.combinatorics.tensor_can import canonical_free
    >>> gens = [[1,0,2,3,5,4], [2,3,0,1,4,5],[0,1,3,2,5,4]]
    >>> gens = [Permutation(h) for h in gens]
    >>> base = [0, 2]
    >>> g = Permutation([2, 1, 0, 3, 4, 5])
    >>> canonical_free(base, gens, g, 4)
    [0, 3, 1, 2, 5, 4]

    Consider the product of Riemann tensors
    `T = R^{a}_{d0}^{d1,d2}*R_{d2,d1}^{d0,b}`
    The order of the indices is [a,b,d0,-d0,d1,-d1,d2,-d2]
    The permutation corresponding to the tensor is
    g = [0,3,4,6,7,5,2,1,8,9]
    Use the slot symmetries to get `T` is the form which is the minimum
    in lexicographic order
    `R^{a}_{d0}^{d1,d2}*R^{b,d0}_{d1,d2}` corresponding to
    `[0, 3, 4, 6, 1, 2, 5, 7, 8, 9]`
    >>> from sympy.combinatorics.tensor_can import riemann_bsgs, tensor_gens
    >>> base, gens = riemann_bsgs
    >>> size, sbase, sgens = tensor_gens(base, gens, [[],[]], 0)
    >>> g = Permutation([0,3,4,6,7,5,2,1,8,9])
    >>> canonical_free(sbase, [Permutation(h) for h in sgens], g, 2)
    [0, 3, 4, 6, 1, 2, 5, 7, 8, 9]
    """
    g = g.array_form
    size = len(g)
    if not base:
        return g[:]

    transversals = get_transversals(base, gens)
    cosets = transversal2coset(size, base, transversals)
    K = [h._array_form for h in gens]
    m = len(base)
    for x in sorted(g[:-2]):
        if x not in base:
            base.append(x)
    lambd = g
    K1 = []
    for i in range(m):
        b = base[i]
        delta = [x[b] for x in cosets[b]]
        delta1 = [lambd[x] for x in delta]
        delta1min = min(delta1)
        k = delta1.index(delta1min)
        p = delta[k]
        for omega in cosets[b]:
            if omega[b] == p:
                break
        lambd = _af_rmul(lambd, omega)
        K = [px for px in K if px[b] == b]
    return lambd