(dummies, sym, b_S, sgens, S_transversals, g)
return TAB[0][-1]
Butler-Portugal algorithm for tensor canonicalization with dummy indices

      dummies
        list of lists of dummy indices,
        one list for each type of index;
        the dummy indices are put in order contravariant, covariant
        [d0, -d0, d1, -d1, ...].

      sym
        list of the symmetries of the index metric for each type.

      possible symmetries of the metrics
              * 0     symmetric
              * 1     antisymmetric
              * None  no symmetry

      b_S
        base of a minimal slot symmetry BSGS.

      sgens
        generators of the slot symmetry BSGS.

      S_transversals
        transversals for the slot BSGS.

      g
        permutation representing the tensor.

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


    A tensor with dummy indices can be represented in a number
    of equivalent ways which typically grows exponentially with
    the number of indices. To be able to establish if two tensors
    with many indices are equal becomes computationally very slow
    in absence of an efficient algorithm.

    The Butler-Portugal algorithm [3] is an efficient algorithm to
    put tensors in canonical form, solving the above problem.

    Portugal observed that a tensor can be represented by a permutation,
    and that the class of tensors equivalent to it under slot and dummy
    symmetries is equivalent to the double coset `D*g*S`
    (Note: in this documentation we use the conventions for multiplication
    of permutations p, q with (p*q)(i) = p[q[i]] which is opposite
    to the one used in the Permutation class)

    Using the algorithm by Butler to find a representative of the
    double coset one can find a canonical form for the tensor.

    To see this correspondence,
    let `g` be a permutation in array form; a tensor with indices `ind`
    (the indices including both the contravariant and the covariant ones)
    can be written as

    `t = T(ind[g[0],..., ind[g[n-1]])`,

    where `n= len(ind)`;
    `g` has size `n + 2`, the last two indices for the sign of the tensor
    (trick introduced in [4]).

    A slot symmetry transformation `s` is a permutation acting on the slots
    `t -> T(ind[(g*s)[0]],..., ind[(g*s)[n-1]])`

    A dummy symmetry transformation acts on `ind`
    `t -> T(ind[(d*g)[0]],..., ind[(d*g)[n-1]])`

    Being interested only in the transformations of the tensor under
    these symmetries, one can represent the tensor by `g`, which transforms
    as

    `g -> d*g*s`, so it belongs to the coset `D*g*S`.

    Let us explain the conventions by an example.

    Given a tensor `T^{d3 d2 d1}{}_{d1 d2 d3}` with the slot symmetries
          `T^{a0 a1 a2 a3 a4 a5} = -T^{a2 a1 a0 a3 a4 a5}`

          `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} -> ...`

    where `[b0, b1, ...] = range(b)` is a base of the symmetric group;
    the strong base `b_S` of S is an ordered sublist of it;
    therefore it is sufficient to compute once the
    strong base generators of S using the Schreier-Sims algorithm;
    the stabilizers of the strong base generators are the
    strong base generators of the stabilizer subgroup.

    `dbase = [p0,p1,...]` is not in general in lexicographic order,
    so that one must recompute the strong base generators each time;
    however this is trivial, there is no need to use the Schreier-Sims
    algorithm for D.

    The algorithm keeps a TAB of elements `(s_i, d_i, h_i)`
    where `h_i = d_i*g*s_i` satisfying `h_i[j] = p_j` for `0 <= j < i`
    starting from `s_0 = id, d_0 = id, h_0 = g`.

    The equations `h_0[0] = p_0, h_1[1] = p_1,...` are solved in this order,
    choosing each time the lowest possible value of p_i

    For `j < i`
    `d_i*g*s_i*S_{b_0,...,b_{i-1}}*b_j = D_{p_0,...,p_{i-1}}*p_j`
    so that for dx in `D_{p_0,...,p_{i-1}}` and sx in
    `S_{base[0],...,base[i-1]}` one has `dx*d_i*g*s_i*sx*b_j = p_j`

    Search for dx, sx such that this equation holds for `j = i`;
    it can be written as `s_i*sx*b_j = J, dx*d_i*g*J = p_j`
    `sx*b_j = s_i**-1*J; sx = trace(s_i**-1, S_{b_0,...,b_{i-1}})`
    `dx**-1*p_j = d_i*g*J; dx = trace(d_i*g*J, D_{p_0,...,p_{i-1}})`

    `s_{i+1} = s_i*trace(s_i**-1*J, S_{b_0,...,b_{i-1}})`
    `d_{i+1} = trace(d_i*g*J, D_{p_0,...,p_{i-1}})**-1*d_i`
    `h_{i+1}*b_i = d_{i+1}*g*s_{i+1}*b_i = p_i`

    `h_n*b_j = p_j` for all j, so that `h_n` is the solution.

    Add the found `(s, d, h)` to TAB1.

    At the end of the iteration sort TAB1 with respect to the `h`;
    if there are two consecutive `h` in TAB1 which differ only for the
    sign, the tensor is zero, so return 0;
    if there are two consecutive `h` which are equal, keep only one.

    Then stabilize the slot generators under `i` and the dummy generators
    under `p_i`.

    Assign `TAB = TAB1` at the end of the iteration step.

    At the end `TAB` contains a unique `(s, d, h)`, since all the slots
    of the tensor `h` have been fixed to have the minimum value according
    to the symmetries. The algorithm returns `h`.

    It is important that the slot BSGS has lexicographic minimal base,
    otherwise there is an `i` which does not belong to the slot base
    for which `p_i` is fixed by the dummy symmetry only, while `i`
    is not invariant from the slot stabilizer, so `p_i` is not in
    general the minimal value.

    This algorithm differs slightly from the original algorithm [3]:
      the canonical form is minimal lexicographically, and
      the BSGS has minimal base under lexicographic order.
      Equal tensors `h` are eliminated from TAB.


    Examples
    ========

    >>> from sympy.combinatorics.permutations import Permutation
    >>> from sympy.combinatorics.perm_groups import PermutationGroup
    >>> from sympy.combinatorics.tensor_can import double_coset_can_rep, get_transversals
    >>> gens = [Permutation(x) for x in [[2,1,0,3,4,5,7,6], [4,1,2,3,0,5,7,6]]]
    >>> base = [0, 2]
    >>> g = Permutation([4,2,0,1,3,5,6,7])
    >>> transversals = get_transversals(base, gens)
    >>> double_coset_can_rep([list(range(6))], [0], base, gens, transversals, g)
    [0, 1, 2, 3, 4, 5, 7, 6]

    >>> g = Permutation([4,1,3,0,5,2,6,7])
    >>> double_coset_can_rep([list(range(6))], [0], base, gens, transversals, g)
    0
Butler-Portugal algorithm for tensor canonicalization with dummy indices
Butler - Portugal algorithm for tensor canonicalization with dummy indices
def double_coset_can_rep(dummies, sym, b_S, sgens, S_transversals, g):
    """
    Butler-Portugal algorithm for tensor canonicalization with dummy indices

      dummies
        list of lists of dummy indices,
        one list for each type of index;
        the dummy indices are put in order contravariant, covariant
        [d0, -d0, d1, -d1, ...].

      sym
        list of the symmetries of the index metric for each type.

      possible symmetries of the metrics
              * 0     symmetric
              * 1     antisymmetric
              * None  no symmetry

      b_S
        base of a minimal slot symmetry BSGS.

      sgens
        generators of the slot symmetry BSGS.

      S_transversals
        transversals for the slot BSGS.

      g
        permutation representing the tensor.

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


    A tensor with dummy indices can be represented in a number
    of equivalent ways which typically grows exponentially with
    the number of indices. To be able to establish if two tensors
    with many indices are equal becomes computationally very slow
    in absence of an efficient algorithm.

    The Butler-Portugal algorithm [3] is an efficient algorithm to
    put tensors in canonical form, solving the above problem.

    Portugal observed that a tensor can be represented by a permutation,
    and that the class of tensors equivalent to it under slot and dummy
    symmetries is equivalent to the double coset `D*g*S`
    (Note: in this documentation we use the conventions for multiplication
    of permutations p, q with (p*q)(i) = p[q[i]] which is opposite
    to the one used in the Permutation class)

    Using the algorithm by Butler to find a representative of the
    double coset one can find a canonical form for the tensor.

    To see this correspondence,
    let `g` be a permutation in array form; a tensor with indices `ind`
    (the indices including both the contravariant and the covariant ones)
    can be written as

    `t = T(ind[g[0],..., ind[g[n-1]])`,

    where `n= len(ind)`;
    `g` has size `n + 2`, the last two indices for the sign of the tensor
    (trick introduced in [4]).

    A slot symmetry transformation `s` is a permutation acting on the slots
    `t -> T(ind[(g*s)[0]],..., ind[(g*s)[n-1]])`

    A dummy symmetry transformation acts on `ind`
    `t -> T(ind[(d*g)[0]],..., ind[(d*g)[n-1]])`

    Being interested only in the transformations of the tensor under
    these symmetries, one can represent the tensor by `g`, which transforms
    as

    `g -> d*g*s`, so it belongs to the coset `D*g*S`.

    Let us explain the conventions by an example.

    Given a tensor `T^{d3 d2 d1}{}_{d1 d2 d3}` with the slot symmetries
          `T^{a0 a1 a2 a3 a4 a5} = -T^{a2 a1 a0 a3 a4 a5}`

          `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} -> ...`

    where `[b0, b1, ...] = range(b)` is a base of the symmetric group;
    the strong base `b_S` of S is an ordered sublist of it;
    therefore it is sufficient to compute once the
    strong base generators of S using the Schreier-Sims algorithm;
    the stabilizers of the strong base generators are the
    strong base generators of the stabilizer subgroup.

    `dbase = [p0,p1,...]` is not in general in lexicographic order,
    so that one must recompute the strong base generators each time;
    however this is trivial, there is no need to use the Schreier-Sims
    algorithm for D.

    The algorithm keeps a TAB of elements `(s_i, d_i, h_i)`
    where `h_i = d_i*g*s_i` satisfying `h_i[j] = p_j` for `0 <= j < i`
    starting from `s_0 = id, d_0 = id, h_0 = g`.

    The equations `h_0[0] = p_0, h_1[1] = p_1,...` are solved in this order,
    choosing each time the lowest possible value of p_i

    For `j < i`
    `d_i*g*s_i*S_{b_0,...,b_{i-1}}*b_j = D_{p_0,...,p_{i-1}}*p_j`
    so that for dx in `D_{p_0,...,p_{i-1}}` and sx in
    `S_{base[0],...,base[i-1]}` one has `dx*d_i*g*s_i*sx*b_j = p_j`

    Search for dx, sx such that this equation holds for `j = i`;
    it can be written as `s_i*sx*b_j = J, dx*d_i*g*J = p_j`
    `sx*b_j = s_i**-1*J; sx = trace(s_i**-1, S_{b_0,...,b_{i-1}})`
    `dx**-1*p_j = d_i*g*J; dx = trace(d_i*g*J, D_{p_0,...,p_{i-1}})`

    `s_{i+1} = s_i*trace(s_i**-1*J, S_{b_0,...,b_{i-1}})`
    `d_{i+1} = trace(d_i*g*J, D_{p_0,...,p_{i-1}})**-1*d_i`
    `h_{i+1}*b_i = d_{i+1}*g*s_{i+1}*b_i = p_i`

    `h_n*b_j = p_j` for all j, so that `h_n` is the solution.

    Add the found `(s, d, h)` to TAB1.

    At the end of the iteration sort TAB1 with respect to the `h`;
    if there are two consecutive `h` in TAB1 which differ only for the
    sign, the tensor is zero, so return 0;
    if there are two consecutive `h` which are equal, keep only one.

    Then stabilize the slot generators under `i` and the dummy generators
    under `p_i`.

    Assign `TAB = TAB1` at the end of the iteration step.

    At the end `TAB` contains a unique `(s, d, h)`, since all the slots
    of the tensor `h` have been fixed to have the minimum value according
    to the symmetries. The algorithm returns `h`.

    It is important that the slot BSGS has lexicographic minimal base,
    otherwise there is an `i` which does not belong to the slot base
    for which `p_i` is fixed by the dummy symmetry only, while `i`
    is not invariant from the slot stabilizer, so `p_i` is not in
    general the minimal value.

    This algorithm differs slightly from the original algorithm [3]:
      the canonical form is minimal lexicographically, and
      the BSGS has minimal base under lexicographic order.
      Equal tensors `h` are eliminated from TAB.


    Examples
    ========

    >>> from sympy.combinatorics.permutations import Permutation
    >>> from sympy.combinatorics.perm_groups import PermutationGroup
    >>> from sympy.combinatorics.tensor_can import double_coset_can_rep, get_transversals
    >>> gens = [Permutation(x) for x in [[2,1,0,3,4,5,7,6], [4,1,2,3,0,5,7,6]]]
    >>> base = [0, 2]
    >>> g = Permutation([4,2,0,1,3,5,6,7])
    >>> transversals = get_transversals(base, gens)
    >>> double_coset_can_rep([list(range(6))], [0], base, gens, transversals, g)
    [0, 1, 2, 3, 4, 5, 7, 6]

    >>> g = Permutation([4,1,3,0,5,2,6,7])
    >>> double_coset_can_rep([list(range(6))], [0], base, gens, transversals, g)
    0
    """
    size = g.size
    g = g.array_form
    num_dummies = size - 2
    indices = list(range(num_dummies))
    all_metrics_with_sym = all([_ is not None for _ in sym])
    num_types = len(sym)
    dumx = dummies[:]
    dumx_flat = []
    for dx in dumx:
        dumx_flat.extend(dx)
    b_S = b_S[:]
    sgensx = [h._array_form for h in sgens]
    if b_S:
        S_transversals = transversal2coset(size, b_S, S_transversals)
    # strong generating set for D
    dsgsx = []
    for i in range(num_types):
        dsgsx.extend(dummy_sgs(dumx[i], sym[i], num_dummies))
    ginv = _af_invert(g)
    idn = list(range(size))
    # TAB = list of entries (s, d, h) where h = _af_rmuln(d,g,s)
    # for short, in the following d*g*s means _af_rmuln(d,g,s)
    TAB = [(idn, idn, g)]
    for i in range(size - 2):
        b = i
        testb = b in b_S and sgensx
        if testb:
            sgensx1 = [_af_new(_) for _ in sgensx]
            deltab = _orbit(size, sgensx1, b)
        else:
            deltab = set([b])
        # p1 = min(IMAGES) = min(Union D_p*h*deltab for h in TAB)
        if all_metrics_with_sym:
            md = _min_dummies(dumx, sym, indices)
        else:
            md = [min(_orbit(size, [_af_new(
                ddx) for ddx in dsgsx], ii)) for ii in range(size - 2)]

        p_i = min([min([md[h[x]] for x in deltab]) for s, d, h in TAB])
        dsgsx1 = [_af_new(_) for _ in dsgsx]
        Dxtrav = _orbit_transversal(size, dsgsx1, p_i, False, af=True) \
            if dsgsx else None
        if Dxtrav:
            Dxtrav = [_af_invert(x) for x in Dxtrav]
        # compute the orbit of p_i
        for ii in range(num_types):
            if p_i in dumx[ii]:
                # the orbit is made by all the indices in dum[ii]
                if sym[ii] is not None:
                    deltap = dumx[ii]
                else:
                    # the orbit is made by all the even indices if p_i
                    # is even, by all the odd indices if p_i is odd
                    p_i_index = dumx[ii].index(p_i) % 2
                    deltap = dumx[ii][p_i_index::2]
                break
        else:
            deltap = [p_i]
        TAB1 = []
        nTAB = len(TAB)
        while TAB:
            s, d, h = TAB.pop()
            if min([md[h[x]] for x in deltab]) != p_i:
                continue
            deltab1 = [x for x in deltab if md[h[x]] == p_i]
            # NEXT = s*deltab1 intersection (d*g)**-1*deltap
            dg = _af_rmul(d, g)
            dginv = _af_invert(dg)
            sdeltab = [s[x] for x in deltab1]
            gdeltap = [dginv[x] for x in deltap]
            NEXT = [x for x in sdeltab if x in gdeltap]
            # d, s satisfy
            # d*g*s*base[i-1] = p_{i-1}; using the stabilizers
            # d*g*s*S_{base[0],...,base[i-1]}*base[i-1] =
            # D_{p_0,...,p_{i-1}}*p_{i-1}
            # so that to find d1, s1 satisfying d1*g*s1*b = p_i
            # one can look for dx in D_{p_0,...,p_{i-1}} and
            # sx in S_{base[0],...,base[i-1]}
            # d1 = dx*d; s1 = s*sx
            # d1*g*s1*b = dx*d*g*s*sx*b = p_i
            for j in NEXT:
                if testb:
                    # solve s1*b = j with s1 = s*sx for some element sx
                    # of the stabilizer of ..., base[i-1]
                    # sx*b = s**-1*j; sx = _trace_S(s, j,...)
                    # s1 = s*trace_S(s**-1*j,...)
                    s1 = _trace_S(s, j, b, S_transversals)
                    if not s1:
                        continue
                    else:
                        s1 = [s[ix] for ix in s1]
                else:
                    s1 = s
                #assert s1[b] == j  # invariant
                # solve d1*g*j = p_i with d1 = dx*d for some element dg
                # of the stabilizer of ..., p_{i-1}
                # dx**-1*p_i = d*g*j; dx**-1 = trace_D(d*g*j,...)
                # d1 = trace_D(d*g*j,...)**-1*d
                # to save an inversion in the inner loop; notice we did
                # Dxtrav = [perm_af_invert(x) for x in Dxtrav] out of the loop
                if Dxtrav:
                    d1 = _trace_D(dg[j], p_i, Dxtrav)
                    if not d1:
                        continue
                else:
                    if p_i != dg[j]:
                        continue
                    d1 = idn
                assert d1[dg[j]] == p_i  # invariant
                d1 = [d1[ix] for ix in d]
                h1 = [d1[g[ix]] for ix in s1]
                #assert h1[b] == p_i  # invariant
                TAB1.append((s1, d1, h1))

        # if TAB contains equal permutations, keep only one of them;
        # if TAB contains equal permutations up to the sign, return 0
        TAB1.sort(key=lambda x: x[-1])
        nTAB1 = len(TAB1)
        prev = [0] * size
        while TAB1:
            s, d, h = TAB1.pop()
            if h[:-2] == prev[:-2]:
                if h[-1] != prev[-1]:
                    return 0
            else:
                TAB.append((s, d, h))
            prev = h

        # stabilize the SGS
        sgensx = [h for h in sgensx if h[b] == b]
        if b in b_S:
            b_S.remove(b)
        _dumx_remove(dumx, dumx_flat, p_i)
        dsgsx = []
        for i in range(num_types):
            dsgsx.extend(dummy_sgs(dumx[i], sym[i], num_dummies))
    return TAB[0][-1]
def double_coset_can_rep ( dummies , sym , b_S , sgens , S_transversals , g ) : size = g . size g = g . array_form num_dummies = size - 2 indices = list ( range ( num_dummies ) ) all_metrics_with_sym = all ( [ _ is not None for _ in sym ] ) num_types = len ( sym ) dumx = dummies [ : ] dumx_flat = [ ] for dx in dumx : dumx_flat . extend ( dx ) b_S = b_S [ : ] sgensx = [ h . _array_form for h in sgens ] if b_S : S_transversals = transversal2coset ( size , b_S , S_transversals ) # strong generating set for D dsgsx = [ ] for i in range ( num_types ) : dsgsx . extend ( dummy_sgs ( dumx [ i ] , sym [ i ] , num_dummies ) ) ginv = _af_invert ( g ) idn = list ( range ( size ) ) # TAB = list of entries (s, d, h) where h = _af_rmuln(d,g,s) # for short, in the following d*g*s means _af_rmuln(d,g,s) TAB = [ ( idn , idn , g ) ] for i in range ( size - 2 ) : b = i testb = b in b_S and sgensx if testb : sgensx1 = [ _af_new ( _ ) for _ in sgensx ] deltab = _orbit ( size , sgensx1 , b ) else : deltab = set ( [ b ] ) # p1 = min(IMAGES) = min(Union D_p*h*deltab for h in TAB) if all_metrics_with_sym : md = _min_dummies ( dumx , sym , indices ) else : md = [ min ( _orbit ( size , [ _af_new ( ddx ) for ddx in dsgsx ] , ii ) ) for ii in range ( size - 2 ) ] p_i = min ( [ min ( [ md [ h [ x ] ] for x in deltab ] ) for s , d , h in TAB ] ) dsgsx1 = [ _af_new ( _ ) for _ in dsgsx ] Dxtrav = _orbit_transversal ( size , dsgsx1 , p_i , False , af = True ) if dsgsx else None if Dxtrav : Dxtrav = [ _af_invert ( x ) for x in Dxtrav ] # compute the orbit of p_i for ii in range ( num_types ) : if p_i in dumx [ ii ] : # the orbit is made by all the indices in dum[ii] if sym [ ii ] is not None : deltap = dumx [ ii ] else : # the orbit is made by all the even indices if p_i # is even, by all the odd indices if p_i is odd p_i_index = dumx [ ii ] . index ( p_i ) % 2 deltap = dumx [ ii ] [ p_i_index : : 2 ] break else : deltap = [ p_i ] TAB1 = [ ] nTAB = len ( TAB ) while TAB : s , d , h = TAB . pop ( ) if min ( [ md [ h [ x ] ] for x in deltab ] ) != p_i : continue deltab1 = [ x for x in deltab if md [ h [ x ] ] == p_i ] # NEXT = s*deltab1 intersection (d*g)**-1*deltap dg = _af_rmul ( d , g ) dginv = _af_invert ( dg ) sdeltab = [ s [ x ] for x in deltab1 ] gdeltap = [ dginv [ x ] for x in deltap ] NEXT = [ x for x in sdeltab if x in gdeltap ] # d, s satisfy # d*g*s*base[i-1] = p_{i-1}; using the stabilizers # d*g*s*S_{base[0],...,base[i-1]}*base[i-1] = # D_{p_0,...,p_{i-1}}*p_{i-1} # so that to find d1, s1 satisfying d1*g*s1*b = p_i # one can look for dx in D_{p_0,...,p_{i-1}} and # sx in S_{base[0],...,base[i-1]} # d1 = dx*d; s1 = s*sx # d1*g*s1*b = dx*d*g*s*sx*b = p_i for j in NEXT : if testb : # solve s1*b = j with s1 = s*sx for some element sx # of the stabilizer of ..., base[i-1] # sx*b = s**-1*j; sx = _trace_S(s, j,...) # s1 = s*trace_S(s**-1*j,...) s1 = _trace_S ( s , j , b , S_transversals ) if not s1 : continue else : s1 = [ s [ ix ] for ix in s1 ] else : s1 = s #assert s1[b] == j  # invariant # solve d1*g*j = p_i with d1 = dx*d for some element dg # of the stabilizer of ..., p_{i-1} # dx**-1*p_i = d*g*j; dx**-1 = trace_D(d*g*j,...) # d1 = trace_D(d*g*j,...)**-1*d # to save an inversion in the inner loop; notice we did # Dxtrav = [perm_af_invert(x) for x in Dxtrav] out of the loop if Dxtrav : d1 = _trace_D ( dg [ j ] , p_i , Dxtrav ) if not d1 : continue else : if p_i != dg [ j ] : continue d1 = idn assert d1 [ dg [ j ] ] == p_i # invariant d1 = [ d1 [ ix ] for ix in d ] h1 = [ d1 [ g [ ix ] ] for ix in s1 ] #assert h1[b] == p_i  # invariant TAB1 . append ( ( s1 , d1 , h1 ) ) # if TAB contains equal permutations, keep only one of them; # if TAB contains equal permutations up to the sign, return 0 TAB1 . sort ( key = lambda x : x [ - 1 ] ) nTAB1 = len ( TAB1 ) prev = [ 0 ] * size while TAB1 : s , d , h = TAB1 . pop ( ) if h [ : - 2 ] == prev [ : - 2 ] : if h [ - 1 ] != prev [ - 1 ] : return 0 else : TAB . append ( ( s , d , h ) ) prev = h # stabilize the SGS sgensx = [ h for h in sgensx if h [ b ] == b ] if b in b_S : b_S . remove ( b ) _dumx_remove ( dumx , dumx_flat , p_i ) dsgsx = [ ] for i in range ( num_types ) : dsgsx . extend ( dummy_sgs ( dumx [ i ] , sym [ i ] , num_dummies ) ) return TAB [ 0 ] [ - 1 ]