xPerm

Fast manipulation of permutation groups

José M. Martín-García
(c) 2003-2020, under GPL

jose@xact.es
http://www.xact.es/
http://groups.google.com/group/xAct

0. Info and copyright

(************************ 0. Info and copyright ***********************)

IntegerPart/@Date[]

xAct`xPerm`$Version={"1.2.3",{2015,8,23}}

$SessionID

0.1. GPL

(* xPerm, a free package for permutation groups in Mathematica *)

(* Copyright (C) 2003-2020 Jose M. Martin-Garcia *)

(* This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation; either version 2of the License,
  or (at your option) any later version.

This program is distributed in the hope that it will be useful, but
 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
 for more details.

You should have received a copy of the GNU General Public License along
 with this program; if not, write to the Free Software Foundation, Inc.,
  59 Temple Place-Suite 330, Boston, MA 02111-1307, USA.
*)

0.2. Info package

(* :Title: xPerm *)

(* :Author: Jose M. Martin-Garcia *)

(* :Summary: Computational group theory for index canonicalization *)

(* :Brief Discussion:
    - Notations for permutations: It is not clear which one is better,
      so that four are encoded: Perm, Images, Cycles and Rules.
    - Permutations act on the right of lists or other permutations.
      This is the choice in GAP and Butler's book.
    - We follow Portugal's notation: permutations map slots to indices.
    - There are two main objectives in this package:
        a) Schreier-Sims algorithm (Butler's book)
        b) choice of coset and double-coset representatives
           (Portugal's papers).
    - This is the Mathematica code of the package. An external C code
      with a subset of capabilities can be linked via MathLink.
    - Mathematica 9 contains the Butler-Portugal algorithm, but xPerm
      does not yet link to it by default.  
*)
  
(* :Context: xAct`xPerm` *)

(* :Package Version: 1.2.3 *)

(* :Copyright: Jose M. Martin-Garcia 2003-2020 *)

(* :History: (see xPerm.History for details)
    - Version 0.1 (May 2003): Schreier-Sims coded from Butler's book.
      Portugal's routines encoded.
    - Version 0.2 (July 2003): Recoding in C, using Images notations,
      and connection through MathLink. Much faster.
    - Version 0.3 (August 2003): Recoded signed permutations. The basic
      mathematical structure is an algebra, and not the {-1,1} x Sn group.
    - Version 0.4 (March 2004): First public release.
    - Version 0.5 (July 2006) 
    - Version 0.6 (August 2007)
    - Version 0.7 (November 2007)
    - Version 1.0 (March 2008): Release for publication.
    - Version 1.1 (September 2010)
    - Version 1.2 (January 2013)
*)

(* :Keywords: permutations, index canonicalization, Butler-Portugal algorithm *)

(* :Source: xPerm.nb *)

(* :Mathematica Version: 5.0 and later *)

(* :Limitations: Only permutation symmetries, not multiterm symmetries *)

0.3. TODO

Link to the new canonicalizer in Mathematica 9, which is generally faster.

1. Begin package

(************************ 1. Begin package ***********************)

{TimeUsed[],SessionTime[],TimeUsed[]/SessionTime[]}

1.1. Begin package

Protect against multiple loading of the package:

With[{xAct`xPerm`Private`xPermSymbols=DeleteCases[Join[Names["xAct`xPerm`*"],Names["xAct`xPerm`Private`*"]],"$Version"|"xAct`xPerm`$Version"]},
Unprotect/@xAct`xPerm`Private`xPermSymbols;
Clear/@xAct`xPerm`Private`xPermSymbols;
]

Decide whether the package is being loaded from other package:

If[Unevaluated[xAct`xCore`Private`$LastPackage]===xAct`xCore`Private`$LastPackage,xAct`xCore`Private`$LastPackage="xAct`xPerm`"];

xAct`xCore`Private`$LastPackage

Temporary fix for the conflicts with Mathematica 8.0:

xAct`xPerm`Cycles;
xAct`xPerm`RightCosetRepresentative;

System`Cycles;
Off[System`Cycles::shdw];

Start the xPerm package. It needs the xCore package. From version 0.6.2 of xCore the package ExpressionManipulation is not automatically read.

BeginPackage["xAct`xPerm`",{"xAct`xCore`"}]

Message:

If[xAct`xCore`Private`$LastPackage=!="Tensors`",
Print[xAct`xCore`Private`bars];
Print["Package xAct`xPerm`  version ",$Version[[1]],", ",$Version[[2]]];
Print["CopyRight (C) 2003-2020, Jose M. Martin-Garcia, under the General Public License."]
];

The Disclaimer of GPL (remove temporarily the General::shdw message to avoid collision with the disclaimer messages of other packages):

Off[General::shdw]
xAct`xPerm`Disclaimer[]:=Print["These are points 11 and 12 of the General Public License:\n\nBECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM `AS IS´ WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.\n\nIN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES."]
On[General::shdw]

If xPerm is not being called from other package then write this GPL short disclaimer:

If[xAct`xCore`Private`$LastPackage==="xAct`xPerm`",
Unset[xAct`xCore`Private`$LastPackage];
Print[xAct`xCore`Private`bars];
Print["These packages come with ABSOLUTELY NO WARRANTY; for details type Disclaimer[]. This is free software, and you are welcome to redistribute it under certain conditions. See the General Public License for details."];
Print[xAct`xCore`Private`bars]]

This is the output of the Disclaimer[] command:

Disclaimer[]

1.2. Usage messages

(* Notations and translation *)
Perm::usage="Perm is the head for permutation lists of the form Perm[{n1, ..., nk}]. Example: Perm[{6, 3, 2, 1, 5, 4}] means that the object in the first place goes to the fourth place and so on. Perm[{}] represents the identity.";
Cycles::usage="Cycles is the head for permutations expressed in disjoint cyclic notation. Singletons are not included. Example: Cycles[{1, 4, 6}, {2, 3}] represents the exchange of points 2 and 3, and the cycle 1->4->6->1. Cycles[] represents the identity.";
Rules::usage="Rules is the head for permutations expressed in rule notation. Singletons are not included. Example: Rules[1->4, 4->6, 6->1, 2->3, 3->2]. Rules[] represents the identity.";
Images::usage="Images is the head for permutations expressed as a list of images. Example: Images[{4, 3, 2, 6, 5, 1}]. Images[{}] represents the identity.";
PermQ::usage="PermQ[ Perm[{n1, ..., nk}] ] yields True if the list {n1, ..., nk} of length k is a rearrangement of the numbers {1, ..., k}. PermQ[ Cycles[cyc1, cyc2, ...] ] gives True if there are no repeated points among cycles cyci. PermQ[ Rules[i1->i2, ...] ] gives True if the rules map a set of integers to itself. PermQ[ Images[{n1, ..., nk}] ] yields True if the list {n1, ..., nk} of length k is a rearrangement of the numbers {1, ..., k}. PermQ[ID] gives True. PermQ gives True on a linear combination of valid permutations. PermQ returns False otherwise.";
ID::usage="ID represents the identity in all notations. ID[perm] gives the identity permutation in the notation used by perm.";
NotationOfPerm::usage="NotationOfPerm[perm] gives the notation of permutation perm. It can be one of {Perm, length}, {Images, length}, Cyles or Rules.";
TranslatePerm::usage="TranslatePerm[perm, notation] translates permutation perm, in any format, to the given notation, which must be one of Perm, {Perm, length}, Cycles, Rules, Images or {Images, length}. TranslatePerm[set, notation] returns the set (with head GenSet, StrongGenSet, Group or Coset) with all its permutations translated to notation.";
RandomPerm::usage="RandomPerm[deg] gives a random permutation in Perm notation. RandomPerm[deg, notation] constructs RandomPerm[deg] and then translates it to the given notation using TranslatePerm[perm, notation].";

(* Operations with permutations *)
PermLength::usage="PermLength[perm] gives the so-called length of permutation perm. If expressed in Perm or Images notation, the length of perm is the length of the list of points. In other cases, it returns PermDeg[perm].";
PermDeg::usage="PermDeg[perm] gives the degree of permutation perm, defined as the largest point moved. The degree of (any form of) the identity is 0. PermDeg[set] gives the largest degree among the permutations of the set (a group, generating set, strong generating set or linear combination).";
PermProduct::usage="PermProduct[perm1, perm2] gives the composition of permutations perm1 and perm2, applied from left to right, that is, the result is equivalent to the application of perm2 after perm1 on a list of points. PermProduct[perm1, perm2, perm3, ...] composes all permutations, assuming that they all use the same notation. PermProduct[perm] gives perm and PermProduct[] gives ID.";
InversePerm::usage="InversePerm[perm] gives the inverse permutation of perm.";
PermPower::usage="PermPower[perm, n] for positive n gives the composition of n copies of perm. For negative n it gives the composition of n copies of InversePerm[perm]. For zero n it gives ID[perm].";
OnPoints::usage="OnPoints[p, perm] yields the image of point p under permutation perm. OnPoints[{p1, ...}, perm] yields the list of images of points pi. Signs of permutations are not considered.";
PermSignature::usage="PermSignature[perm] gives the signature of perm: +1 if perm can be written as an even number of transpositions, -1 if perm can be written as an odd number of transpositions.";

(* Stabilization *)
StablePoints::usage="StablePoints[perm, n] gives the points in {1, ..., n} that are fixed by permutation perm. StablePoints[GS, n] gives the points in {1, ..., n} that are fixed by all permutations in GS.";
NonStablePoints::usage="NonStablePoints[list, GS] enlarges list until none of the permutations in GS fixes all points in list. It tries to give the shortest result possible, not necessarily ordered.";
Stabilizer::usage="Stabilizer[{p1, ...}, GS] gives the subset of permutations in GS that fix all points pi. Stabilizer[{p1, ...}, SGS] returns a strong generating set for the stabilizer subgroup of the group SGS fixing all points pi. In this latter case a third argument can be given with the length of the internal Schreier vectors.";
SetStabilizer::usage="SetStabilizer[{p1, ...}, GS] gives the subset of permutations in GS that keep the points pi inside the list. SetStabilizer[{p1, ...}, SGS] returns a strong generating set for the setwise stabilizer of the list of points pi under the group generated by the strong generating set SGS.";

(* Permuting lists *)
PermutationFromTo::usage="PermutationFromTo[list1, list2] gives the permutation (in Images notation) that converts list1 into list2. Points in list1 or list2 can be repeated.";
PermuteList::usage="PermuteList[list, perm] permutes the points of list according to the permutation perm.";

(* Sorting *)
SortB::usage="SortB[{e1, e2, ...}, B] sorts the elements ei according to the canonical order defined by list B. Elements not in B are sorted last, using Sort.";
MinB::usage="MinB[list, B] gives the mininum element of list with respect to base B. If none of the elements of B is in list then the minimum of list with respect to the canonical order is returned.";
PermOrderedQ::usage="PermOrderedQ[h[perm1, perm2]] gives True if perm1 maps integers to smaller points than perm2. It gives False if perm1 maps integers to larger points than perm2. It gives Null if perm1 and perm2 maps integers to the same points. PermOrderedQ[h[perm1, perm2], B] gives True if perm1 maps the points of base B to smaller (according to B) images than perm2. It gives False if perm1 maps the points of B to larger (according to B) images than perm2. It gives Null if perm1 and perm2 map the points of B to the same images.";
PermSort::usage="PermSort[{perm1, perm2, ...}] sorts permutations permi according to PermOrderedQ. PermSort[{perm1, perm2, ...}, B] sorts permutations permi according to PermOrderedQ using base B.";
PermEqual::usage="PermEqual[perm1, perm2] gives True if both permutations map all points to the same images, and False otherwise.";
PermLess::usage="PermLess[perm1, perm2] gives True if perm1 maps sorted integers to smaller points than perm2, and False otherwise.";
PermGreater::usage="PermGreater[perm1, perm2] gives True if perm1 maps sorted integers to larger points than perm2, and False otherwise.";
PermLessEqual::usage="PermLessEqual[perm1, perm2] gives True if perm1 maps sorted integers to smaller or equal points than perm2, and False otherwise.";
PermGreaterEqual::usage="PermGreaterEqual[perm1, perm2] gives True if perm1 maps sorted integers to larger of equal points than perm2, and False otherwise.";

(* Structures *)
GenSet::usage="GenSet is the head for a generating set.";
Group::usage="Group is the head for a group of permutations.";
Coset::usage="Coset is the head for a coset of permutations.";
Schreier::usage="Schreier is the head for a Schreier vector. The format is Schreier[orbits, nu, w], where nu is the actual Schreier vector and w is the associated vector of backward pointers.";
StrongGenSet::usage="StrongGenSet is the head for a strong generating set. The format is StrongGenSet[base, GS].";
Dimino::usage="Dimino[GenSet[perm1, ...]] gives the group generated by permutations permi, using Dimino's algorithm.";

(* Strong generating sets *)
Orbit::usage="Orbit[p, GS] gives the orbit of point p under the group generated by GS. The orbit is given as a (not necessarily ordered) list having p as first point.";
Orbits::usage="Orbits[GS, len] gives a list with all orbits of the group generated by permutations in GS, a partition of Range[len]. Orbits[GS, len, list] makes points in list to be the first points of the orbits. Orbits[GS] and Orbits[GS, list] compute len from GS using PermLength. GS can be a generating set (head GenSet) or a strong generating set (head StrongGenSet).";
SchreierOrbit::usage="SchreierOrbit[p, GS, len] gives the orbit of point p under the group generated by GS, and its Schreier vector (of lenght len). SchreierOrbit[p, SGS, len] gives the same output, but the the first point of the orbit is now the minimum with respect to the base of the SGS. SchreierOrbit[p, gs] computes the length of the vectors from gs using PermLength. See notes for Schreier.";
SchreierOrbits::usage="SchreierOrbits[GS, len] gives a compound Schreier vector (length len) of the orbits of the group generated by GS. SchreierOrbits[GS, len, list] makes points in list to be the first points of the orbits. SchreierOrbits[SGS, len] gives the same output, but the orbits now start with the points of the base of the SGS. SchreierOrbits[gs] computes len from gs using PermLength. See notes for Schreier.";
TraceSchreier::usage="TraceSchreier[p, Schreier[orbits, nu, w]] finds a permutation such that the image of the first point of the orbit of p is p.";
StabilizerChain::usage="StabilizerChain[SGS] gives a list of strong generating sets corresponding to the chain of stabilizers of the group generated by SGS. The first one is SGS and the last one is the generating set of Group[ID].";
OrderOfGroup::usage="OrderOfGroup[SGS] gives the order of the group generated by SGS.";
PermMemberQ::usage="PermMemberQ[perm, SGS] gives True if permutation perm belongs to the group generated by SGS, and False otherwise.";
PermWord::usage="PermWord[perm, SGS] decompose the permutation perm as {pk, ..., p0} where each of the permutations pi belong to the stabilizer subgroup of the first i points of the base of SGS. If pk == ID then perm belongs to the group and not otherwise. There is a third argument, used internally for iteration.";
FromBaseImage::usage="FromBaseImage[list, SGS] gives the permutation in the group described by SGS such that the images of the points in the base are the points in list (assumed to have the length of the base of the SGS). If no permutation in the group gives the list of points then an error message is thrown. A third argument can be given with the length of the internal Schreier vectors used.";
AllBaseImages::usage="AllBaseImages[SGS] gives a complete list of rules baseimages -> perm, where baseimages is a list of images of the base of SGS and perm is its corresponding permutation, as would be given by FromBaseImage[baseimages, SGS]. A second argument can be given with the length of the Schreier vectors used.";
Search::usage="Search[SGS, P, s, SGSK] returns the s-th stabilizer in the stabilizer-chain of the subgroup K of permutations obeying the property P in the group G (described by the strong generating set SGS). It is assumed P[perm] returns True or False on any permutation of the group G. The fourth argument is a strong generating set for a subgroup of K, possibly a deeper stabilizer in its stabilizer-chain.";
BaseChange::usage="BaseChange[SGS, newbase] gives a new strong generating set for the group described by SGS, having newbase as the first elements of its base. A third argument can be given specifying the length of the internal Schreier vectors used.";
DeleteRedundantGenerators::usage="DeleteRedundantGenerators[SGS] returns an equivalent strong generating set removing redundant generators from the generating set.";
SchreierSims::usage="SchreierSims[initbase, GS, deg] generates a strong generating set for the group generated by GS (permutations of degree deg), using list initbase as the first points for the base. The final SGS is not reduced in general.";
UseRules::usage="UseRules is an option for SchreierSims giving a set of rules replacing the permutations by strings. By default it is {}.";
xPermVerbose::usage="xPermVerbose is an option for SchreierSims, RightCosetRepresentative, DoubleCosetRepresentative and CanonicalPerm. xPermVerbose->True gives lots of information about the intermediate status of the process.";

(* Canonicalization *)
RightCosetRepresentative::usage="RightCosetRepresentative[perm, n, SGS] for an unsigned permutation perm of degree n and a SGS for group S gives a canonical representative of the right coset S.perm of perm with respect to subgroup S of the symmetric group Sn. The criterium is the minimization of images of points of the base of SGS under the elements of the coset, following the order given by the base. A fourth argument can be used to give additional priority to some points (the free slots in the tensorial context)..";
DummySet::usage="DummySet is the head of expressions DummySet[manifold, {{d1u, d1d}, {d2u, d2d}, ...}, metricsym] denoting a set of pairs of (names of) dummies {diu, did} on manifold, whose metric has symmetry metricsym (an integer with value 1 if the metric is symmetric, -1 if the metric is antisymmetric, or 0 if there is no metric).";
RepeatedSet::usage="RepeatedSet[{i1, i2, ...}] represents a list of names of repeated indices in the canonical configuration.";
JoinSGS::usage="JoinSGS[StrongGenSet[base1, GS1], StrongGenSet[base2, GS2]] gives a strong generating set having base Join[base1, base2] (with base1 and base2 assumed to be disjoint) and generating set Union[GS1, GS2] (with GS1 and GS2 assumed to move disjoint sets of points).";
SGSOfDummySet::usage="SGSOfDummySet[ DummySet[...] ] gives a Strong Generating Set for the group of permutations associated to the given DummySet. There are always permutations coming from the exchange of dummies. There are permutations coming from the exchange of up/down indices in a pair if there is a metric for those indices. The GS is given in Cycles notation. See notes for DummySet.";
DoubleCosetRepresentative::usage="DoubleCosetRepresentative[perm, n, SGS, dummysets] returns a canonical representative of the double coset S.perm.D, where S is generated by the strong generating set SGS and D is the group of symmetries of the dummysets. The criterium is the sequential sorting of slots of permutation perm (from first to last), choosing the least index not yet used consistent with the symmetries S and D. The indices are chosen assuming that indices are numbered according to their canonical order.";
CanonicalPerm::usage="CanonicalPerm[perm, n, GS, {f1, f2, ...}, dummysets] gives a canonical representative of the double coset S.perm.D, where S is generated by GS and D is the group of symmetries of the dummysets of dummies. The algorithm RightCosetRepresentative is first applied to canonicalize free indices f1, f2, ..., and then algorithm DoubleCosetRepresentative is applied on the result of the latter.";
OrderedBase::usage="OrderedBase is an option for CanonicalPerm specifying that the base must be filled with the missing integers, in normal order.";
BaseChangeCheck::usage="BaseChangeCheck is a Boolean option for RightCosetRepresentative that forces checking the need of a change of base to stabilize the group of symmetry slots. The default is False.";
TimeVerbose::usage="TimeVerbose is an option for CanonicalPerm. TimeVerbose->True outputs the time spent in the computation of the SGS, in the computation of the free-indices representative and in the computation of the dummy-indices representative.";
MathLink::usage="MathLink is a Boolean option for CanonicalPerm, SchreierSims and Orbit, indicating whether the external C executable xperm must be used to speed up computations. Its default value is given by the global variable $xpermQ for CanonicalPerm, and it is False for the other two functions.";

(* MathLink connection *)
$xpermQ::usage="$xpermQ is a Boolean global variable containing whether the connection to the external C executable xperm has been possible or not.";
$xpermExecutable::usage="$xpermExecutable is a global variable containing the name of the xperm executable to which we have connected. If the connection was not possible then it is not assigned.";
$xpermLink::usage="$xpermLink is a global variable containing the link identification of the connection to the external excutable xperm. If the connection was not possible then it is not assigned.";

(* Predefined strong generating sets *)
If[System`$VersionNumber<8.5,
Symmetric::usage="Symmetric[{p1, p2, ...}] represents the symmetry of a tensor that is symmetric in the slots pi. Symmetric[{p1, p2, ...}, notation] returns a strong generating set for the symmetric group on the set of points pi, using the indicated notation for signed permutations.";
Antisymmetric::usage="Antisymmetric[{p1, p2, ...}] represents the symmetry of a tensor that is antisymmetric in the slots pi. Antisymmetric[{p1, p2, ...}, notation] returns a strong generating set for the antisymmetric group on the set of points pi, using the indicated notation for signed permutations.";,
Symmetric::usage=Symmetric::usage<>"\n\nxAct extension:\nSymmetric[{p1, p2, ...}, notation] returns a strong generating set for the symmetric group on the set of points pi, using the indicated notation for signed permutations.";
Antisymmetric::usage=Antisymmetric::usage<>"\n\nxAct extension:\nAntisymmetric[{p1, p2, ...}, notation] returns a strong generating set for the antisymmetric group on the set of points pi, using the indicated notation for signed permutations.";
];
PairSymmetric::usage="PairSymmetric[{{p1a, p1b}, {p2a, p2b}, ...}, sym1, sym2, notation] returns a Strong Generating Set for the symmetric (if sym1=1) or antisymmetric (if sym1=-1) group of permutations of the pairs {pia, pib}. The switch sym is an integer: 1 adds permutations Cycles[{pia, pib}]; -1 adds permutations -Cycles[{pia, pib}]; other values do nothing. The result is given in the using the indicated notation (Cycles by default).";
RiemannSymmetric::usage="RiemannSymmetric[{p1, p2, p3, p4}, notation] gives a strong generating set implementing the symmetries of the Riemann tensor R_{p1,p2,p3,p4}, using the indicated notation for signed permutations (Cycles by default).";
RiemannSymmetry::usage="RiemannSymmetry = RiemannSymmetric. Kept for backward compatibility.";

(* Transversals *)
RightTransversal::usage="RightTransversal[SGS, deg] returns a sorted list of permutations containing the canonical representative of each right coset of SGS in the symmetric group of degree deg. The choice of representative is based in the order induced on the permutations by the base of the SGS. RightTransversal[GS, deg] returns the same thing by first using SchreierSims.";
LeftTransversal::usage="LeftTransversal[SGS, deg] returns a sorted list of permutations containing the canonical representative of each left coset of SGS in the symmetric group of degree deg. The choice of representative is based in the order induced on the **inverse** of the permutations by the base of the SGS. Efectively the transversal is simply the inverse of the corresponding right transversal. LeftTransversal[GS, deg] returns the same thing by first using SchreierSims.";
DoubleTransversal::usage="DoubleTransversal[SGS, dummysets] returns a sorted list of permutations containing the canonical representative of each double coset of SGS and the D-group of the dummysets in the symmetric group of degree deg.";

1.3. Begin Private

Begin["`Private`"]

$xPermNames=Names["xAct`xPerm`*"];

Length[$xPermNames]

There are 76 public reserved words:

?xAct`xPerm`*

2. Permutation Defs. Notation dependent

(************** 2. Permutations. Notation dependent **************)

{TimeUsed[],SessionTime[],TimeUsed[]/SessionTime[]}

2.0. Comments. Relation with Mathematica’s GroupTheory code

If[$ReadingVerbose,Print["Reading section 2: Permutations. Notation dependent"],Null,Null]

This section defines all those actions on permutations whose structure is different on different notations for the permutations.

2.0.1. Fake code in Mathematica 8 for pre8 versions

If[System`$VersionNumber<8.,

PermutationListQ[list_]:=SameQ[Sort[list],Range[Length[list]]];

PermutationCyclesQ[System`Cycles[{cycs___}]]:=Union[cycs]===Sort[Join[cycs]];

lastFalse[{___,{pos_}}]:=pos;
lastFalse[{}]:=0;
PermutationMax[list_List]:=lastFalse[Position[Inner[SameQ,list,Sort[list],List],False]];

InversePermutation[list_List]:=Ordering[list];
InversePermutation[System`Cycles[cyclist_List]]:=System`Cycles[Reverse/@cyclist];
]

2.0.2. Reuse code from System`. Beware of the opposite conventions in Mma8 and Mma9

Definitions are given so that they follow the conventions in Mathematica 9.

SetAttributes[IfMathematica789,HoldAll];
IfMathematica789[expr7_,expr8_,expr9_]:=Which[
System`$VersionNumber<7.5,expr7,
System`$VersionNumber<8.5,expr8,
True,expr9
];

IfMathematica789[
(* Mathematica 7*)
addsingletons[cycs_,length_]:=Join[cycs,Partition[Complement[Range@length,Flatten@cycs],1]];PermutationFromCycles[cycs_List]:=Last/@Sort@Transpose[Flatten/@{RotateRight/@cycs,cycs}];
PermList[System`Cycles[cycslist_List],length_]:=PermutationFromCycles[addsingletons[cycslist,length]],
(* Mathematica 8 *)
PermList[System`Cycles[cycslist_List], length_]:=PermutationList[System`Cycles[Reverse/@cycslist],length],
(* Mathematica 9 *)
PermList[perm_System`Cycles,length_]:=PermutationList[perm,length]
]

IfMathematica789[
(* Mathematica 7*)
PermList[list_List,length_]:=Take[list~Join~Range[Length[list]+1,length],length],
(* Mathematica 8 *)
PermList[list_List, length_]:=PermutationList[list,length],
(* Mathematica 9 *)
PermList[list_List,length_]:=PermutationList[list,length]
]

IfMathematica789[
(* Mathematica 7 *)
ToCycles[list_List]:=Take[#,Position[Rest@#,First@#][[1,1]]]&/@Fold[If[MemberQ[Flatten@#1,#2],#1,Append[#1,NestList[list[[#1]]&,#2,Length@list]]]&,{},list];
PermCycles[list_List]:=System`Cycles[Reverse/@ToCycles[Ordering@list]],
(* Mathematica 8 *)
PermCycles[list_List]:=Map[Reverse,PermutationCycles[list],{2}],
(* Mathematica 9 *)
PermCycles[list_List]:=PermutationCycles[list]
]

Note that this is only defined for System`Cycles, for which conventions agree in Mma8 and Mma9:

IfMathematica789[
(* Mathematica 7 *)
PermutationProduct[System`Cycles[cycslist1_],System`Cycles[cycslist2_]]:=System`Cycles[List@@TranslatePerm[PermProduct[TranslatePerm[Cycles@@cycslist1,Rules],TranslatePerm[Cycles@@cycslist2,Rules]],Cycles]],
(* Mathematica 8 *)
Null,
(* Mathematica 9 *)
Null
]

2.0.3. Tests

This could be moved to xCore`.

Clear[test,testcounter,testresult];
SetAttributes[test,HoldAll];
testcounter[symbol_]:=(testcounter[symbol]=0);
test[symbol_,expr_,result_,messageQ_:False]:=
{symbol,1+testcounter[symbol]++,testresult[symbol,testcounter[symbol]]=Check[SameQ[expr,result],messageQ]};

alltests[]:=Last/@DownValues[testresult];
alltests[symbol_]:=
Cases[DownValues[testresult],_[_[testresult[symbol,_]],result_]:>result]

2.0.4. Notation conversion to/from Mathematica’s GroupTheory code

From Mathematica to xPerm notation. Note that SchreierSims can give a nonordered basis here. We need to find a way to avoid this call of SchreierSims.

MathToxPermSym[sym:(_Symmetric|_Antisymmetric)]:=SchreierSims[{},sym];
MathToxPermSym[gs:{___List}]:=SchreierSims[{},GenSet@@(MathToxPermGen/@gs)];

MathToxPermGen[{System`Cycles[{cycs___}],sign_}]:=sign Cycles[cycs];

From xPerm to Mathematica notation:

xPermToMathSym[sym:(_Symmetric|_Antisymmetric)]:=sym;
xPermToMathSym[StrongGenSet[base_List,GS_GenSet]]:={xPermToMathGen/@List@@GS,base};

xPermToMathGen[sign_. Cycles[cycs___]]:={System`Cycles[{cycs}],sign};

2.1. Notations for permutations. PermQ

2.1.1. Definitions

This is the function in charge of certifying that something is a permutation:

?PermQ

Identity:

?ID

PermQ[ID]:=True;

Perm notation: A permutation is any of the n! rearrangements of the list of numbers {1,2,...,n}

?Perm

PermQ[Perm[list:{___Integer}]]:=PermutationListQ[list];

Images notation: A permutation is any of the n! rearrangements of the list of numbers {1,2,...,n}

?Images

PermQ[Images[list:{___Integer}]]:=PermutationListQ[list];

Cycles notation: A permutation is any set of disjoint cycles (lists of different points)

?Cycles

PermQ[Cycles[cycs:{___Integer?Positive}...]]:=PermutationCyclesQ[System`Cycles[{cycs}]];

Rules notation: A permutation is any collection of rules mapping a set of integers to itself

?Rules

PermQ[Rules[rules:Rule[_Integer?Positive,_Integer?Positive]..]]:=Apply[Union[#1]===Sort[Join[#2]]&,Transpose[{rules}/.Rule->List]
];PermQ[Rules[]]:=True;

Action on a linear combination. This is not an evident definition; here we have in mind a future generalization to the algebra of permutations:

PermQ[expr_Plus]:=Apply[And,PermQ/@List@@expr];
PermQ[x_?NumericQ perm_]:=PermQ[perm];

Any other case:

PermQ[_]:=False;

Protect:

SetNumberOfArguments[PermQ,1];
Protect[PermQ];

2.1.2. Examples

In the following we do not check that arguments of functions are actually permutations, because that would slow down the computations. The function PermQ is given to test permutations in case of doubt.

Examples:

test[PermQ,
PermQ[-Perm[{1,2,3}]+2Perm[{3,4,5,2,1}]],
True
]

test[PermQ,
PermQ[-Perm[{1,2,3}]Perm[{3,4,5,2,1}]],
False
]

Perm[{}] represents the identity permutation.

test[PermQ,
PermQ[Perm[{}]],
True
]

test[PermQ,
PermQ[Perm[]],
False
]

test[PermQ,
PermQ[Perm[{3,2,1}]],
True
]

test[PermQ,
PermQ[Perm[{1,2,1}]],
False
]

Images[{}] represents the identity permutation.

test[PermQ,
PermQ[Images[{}]],
True
]

test[PermQ,
PermQ[Images[]],
False
]

test[PermQ,
PermQ[Images[{3,2,1}]],
True
]

test[PermQ,
PermQ[Images[{1,2,1}]],
False
]

Cycles:

test[PermQ,
PermQ[Cycles[]],
True
]

test[PermQ,
PermQ[Cycles[{1,3},{2,5}]],
True
]

test[PermQ,
PermQ[Cycles[{1,3},{2,5,0}]],
False
]

test[PermQ,
PermQ[Cycles[{1,3.},{2,5}]],
False
]

test[PermQ,
PermQ[Cycles[{1,3},{1,5}]],
False,
True
]

test[PermQ,
PermQ[Cycles[{1,3,1},{2,5}]],
False,
True
]

Rules:

test[PermQ,
PermQ[Rules[]],
True
]

test[PermQ,
PermQ[Rules[1->3,3->5,5->1]],
True
]

test[PermQ,
PermQ[Rules[1->2,1->3,2->3]],
False
]

And@@alltests[PermQ]

2.1.3. Remove singletons in Cycles and Rules notations

We automatically remove singletons and empty cycles from the cyclic structure (I've tried to replace these definitions with a function DeleteSingletons to be used at the right places, but it is less efficient!)

Cycles[cycs1___,{},cycs2___]:=Cycles[cycs1,cycs2];
Cycles[cycs1___,{_},cycs2___]:=Cycles[cycs1,cycs2];

And the same for rules:

Rules[rules1___,HoldPattern[x_->x_],rules2___]:=Rules[rules1,rules2];

Examples:

Cycles[{},{1},{2,3}]

Rules[1->2,2->2,2->3,3->4]

Note that rules are not sorted:

Rules[2->3,1->2,3->4]

Protect heads:

SetNumberOfArguments[Perm,1];
SetNumberOfArguments[Images,1];
Protect[Perm,Cycles,Rules,Images];

It is important to note that a given permutation has a unique representation in Perm or Images notation (for a given list-length), but not in Cycles or Rules notation because the cycles and rules can be reordered in different ways. It would be simple to solve the problem in Rules notation by introducing the Orderless attribute, but then the computations would slow down. Even worse, it is not so simple for Cycles and actually we do not do it for neither of them. See section 3.1 for tests of equality of permutations.

2.1.4. More on Cycles notation. TOCOMPLETE

Note that the opposite convention for multiplication would have made this nicer :

PermProduct[Cycles[{5,6,7,8,9}],Cycles[{1,2,3,4,5}]]

We have checked that this is the result in GAP as well :

PermProduct[Cycles[{1,2,3,4,5}],Cycles[{5,6,7,8,9}]]

Ensure that cycles are properly sorted:

SortCycles[perm_Cycles]:=Sort[SortCycle/@perm];
SortCycle[cyc_List]:=Nest[RotateLeft,cyc,Position[cyc,Min[cyc]][[1,1]]-1];
SortCycles[expr_]:=expr/.perm_Cycles:>SortCycles[perm];

-3Cycles[{4,3,2,5},{6,1,7,8}]//SortCycles

Standard formatting for cyclic permutation. Note the exception of the identity.

Unprotect[Cycles];
MakeBoxes[Cycles[],StandardForm]:=xAct`xTensor`Private`interpretbox[Cycles[],"id"];
formatcycle[cyc_List]:=RowBox@Join[{"("},Riffle[MakeBoxes/@cyc,","],{")"}];
MakeBoxes[-Cycles[cycs__List],StandardForm]:=xAct`xTensor`Private`interpretbox[-Cycles[cycs],RowBox[Prepend[formatcycle/@{cycs},"-"]]];
MakeBoxes[coeff_?NumberQ Cycles[cycs__List],StandardForm]:=xAct`xTensor`Private`interpretbox[coeff Cycles[cycs],RowBox[Prepend[formatcycle/@{cycs},MakeBoxes[coeff,StandardForm]]]];
MakeBoxes[Cycles[cycs__List],StandardForm]:=xAct`xTensor`Private`interpretbox[Cycles[cycs],RowBox[formatcycle/@{cycs}]];
Protect[Cycles];

Examples :

-Cycles[]

-Cycles[{1,2},{3,4}]

-3Cycles[{1,2},{3,4}]

-3/2Cycles[{1,2},{3,4}]

2.2. Degree and length of permutations. PermDeg, PermLength, NotationOfPerm

2.2.1. PermDeg and PermLength

The length of a permutation is a notation-dependent concept. This is mostly used internally to avoid recomputing the degree of a permutation many times. Users should use PermDeg, which gives the true degree of a permutation, i.e., the largest point moved by the permutation (a notation-independent concept).

?PermLength

We need definitions on the two special cases:

PermLength[Perm[list_List]]:=Length[list];
PermLength[Images[list_List]]:=Length[list];

And definitions on sets of permutations. We take into account the base to compute the perm-length of a strong generating set.

PermLength[group_Group]:=Max[0,PermLength/@List@@group]
PermLength[GS_GenSet]:=Max[0,PermLength/@List@@GS];
PermLength[StrongGenSet[base_,GS_]]:=Max[base,PermLength[GS]];
PermLength[Symmetric[list_]]:=Max[0,list];
PermLength[Antisymmetric[list_]]:=Max[0,list];

In other cases change to PermDeg:

PermLength[x_]:=PermDeg[x];
SetNumberOfArguments[PermLength,1];
Protect[PermLength];

We define PermDeg here, because even though the concept is notation-independent, the internal definitions are different for different notations, to make it faster:

?PermDeg

PermDeg[Perm[list_List]]:=PermutationMax[list];
PermDeg[Images[list_List]]:=PermutationMax[list];
PermDeg[Cycles[cycs___List]]:=Max[0,cycs];
PermDeg[Rules[rules___Rule]]:=Max[0,Apply[List,{rules},{1}]];
PermDeg[ID]:=0;

We do not take into account the base to compute the degree:

PermDeg[Group[perms___]]:=Max[0,PermDeg/@{perms}];
PermDeg[GenSet[perms___]]:=Max[0,PermDeg/@{perms}];
PermDeg[StrongGenSet[_List,GenSet[perms___]]]:=Max[0,PermDeg/@{perms}];
PermDeg[Symmetric[list_]]:=Max[0,list];
PermDeg[Antisymmetric[list_]]:=Max[0,list];

Again, this is just a tentative definition on linear combinations of permutations:

PermDeg[expr_Plus]:=Max[PermDeg/@List@@expr];
PermDeg[x_?NumericQ perm_]:=PermDeg[perm];

PermDeg[x_]:=Throw[Message[PermDeg::undef,"degree",x]];
SetNumberOfArguments[PermDeg,1];
Protect[PermDeg];

2.2.2. Examples

Note the difference between these two cases:

test[PermDeg,
PermLength[Perm[{1,2,3,4}]],
4
]

test[PermDeg,
PermDeg[Perm[{1,2,3,4}]],
0
]

or these:

test[PermDeg,
PermLength[Images[{2,3,1,4}]],
4
]

test[PermDeg,
PermDeg[Images[{2,3,1,4}]],
3
]

Special cases:

test[PermDeg,
PermDeg[Perm[{}]],
0
]

test[PermDeg,
Catch@PermDeg[Perm[]],
Null,
True
]

test[PermDeg,
PermDeg[Rules[]],
0
]

test[PermDeg,
PermDeg[Images[{}]],
0
]

Cycles and Rules notations:

test[PermDeg,
PermDeg[Cycles[{1,2},{4,5}]],
5
]

test[PermDeg,
PermDeg[Cycles[]],
0
]

test[PermDeg,
PermDeg[Rules[1->2,2->5,5->1]],
5
]

And@@alltests[PermDeg]

2.2.3. NotationOfPerm

This simple function gives the notation in which a particular permutation is written. Note that it only acts on single (possibly signed) permutations.

?NotationOfPerm

Note that for a generating set the NotationOfPerm is defined as the notation of the first permutation of the set:

NotationOfPerm[-g_]:=NotationOfPerm[g];
NotationOfPerm[g:Perm[_List]]:={Perm,PermLength[g]};
NotationOfPerm[g:Images[_List]]:={Images,PermLength[g]};
NotationOfPerm[Cycles[___List]]:=Cycles;
NotationOfPerm[Rules[___Rule]]:=Rules;
NotationOfPerm[GenSet[]]:=Cycles;
NotationOfPerm[GenSet[perm_,___]]:=NotationOfPerm[perm];
NotationOfPerm[StrongGenSet[_,GS_]]:=NotationOfPerm[GS];
(* Arbitrary decision *)
NotationOfPerm[Symmetric[list_]]:=Cycles;
NotationOfPerm[Antisymmetric[list_]]:=Cycles;

NotationOfPerm[x_]:=Throw[Message[NotationOfPerm::undef,"notation",x]];
SetNumberOfArguments[NotationOfPerm,1];
Protect[NotationOfPerm];

Examples:

test[NotationOfPerm,
NotationOfPerm[Perm[{3,2,1}]],
{Perm,3}
]

test[NotationOfPerm,
NotationOfPerm[Images[{3,2,1,4}]],
{Images,4}
]

test[NotationOfPerm,
NotationOfPerm[Cycles[{5,4}]],
Cycles
]

test[NotationOfPerm,
NotationOfPerm[Rules[1->100,100->1]],
Rules
]

Catch@NotationOfPerm[Hello]

And@@alltests[NotationOfPerm]

2.3. Identity. ID

2.3.1. Definitions

?ID

Identity permutation in the notation given. Note that the length of the permutation is lost in Cycles and Rules notation:

ID[Perm[list_List]]:=Perm[Sort@list];
ID[Images[list_List]]:=Images[Sort@list];
ID[Cycles[___List]]:=Cycles[];
ID[Rules[___Rule]]:=Rules[];
ID[ID]:=ID;

It is not clear how to define ID on a sum, and so we do not do that.

ID[x_?NumericQ perm_]:=ID[perm];

Protection:

ID[x_]:=Throw[Message[ID::undef,"identity for notation",x]];
SetNumberOfArguments[ID,1];
Protect[ID];

2.3.2. Examples

test[ID,
ID[Perm[{4,3,2,1}]],
Perm[{1,2,3,4}]
]

test[ID,
ID[Perm[{}]],
Perm[{}]
]

test[ID,
ID[Cycles[{2,3},{1,5}]],
Cycles[]
]

test[ID,
ID[-Cycles[{2,3},{1,5}]],
Cycles[]
]

Catch@Perm[]

And@@alltests[ID]

2.4. Notation translations. TranslatePerm

2.4.1. General definitions

This function translates any valid permutation in xPerm to any valid notation.

?TranslatePerm

Thread translate over lists and sets of permutations:

TranslatePerm[set:(_List|_GenSet|_Group|_Coset),notation_]:=TranslatePerm[#,notation]&/@set;
TranslatePerm[StrongGenSet[base_List,GS_GenSet],notation_]:=StrongGenSet[base,TranslatePerm[GS,notation]];
TranslatePerm[Symmetric[list_],notation_]:=Symmetric[list];
TranslatePerm[Antisymmetric[list_],notation_]:=Antisymmetric[list];

Linear combinations:

TranslatePerm[expr_Plus,notation_]:=TranslatePerm[#,notation]&/@expr;
TranslatePerm[x_?NumericQ perm_,notation_]:=x TranslatePerm[perm,notation];

Translation to Perm and Images notations is given in terms of {Perm, length} and {Images, length}:

TranslatePerm[perm_,Perm]:=TranslatePerm[perm,{Perm,PermLength@perm}];
TranslatePerm[perm_,Images]:=TranslatePerm[perm,{Images,PermLength@perm}];

2.4.2. ID to ...

TranslatePerm[ID,{Perm,length_Integer}]:=Perm[Range[length]];
TranslatePerm[ID,{Images,length_Integer}]:=Images[Range[length]];
TranslatePerm[ID,Cycles]:=Cycles[];
TranslatePerm[ID,Rules]:=Rules[];

test[TranslatePerm,
TranslatePerm[ID,Perm],
Perm[{}]
]

test[TranslatePerm,
TranslatePerm[-ID,Perm],
-Perm[{}]
]

test[TranslatePerm,
TranslatePerm[ID,Images],
Images[{}]
]

test[TranslatePerm,
TranslatePerm[-ID,Images],
-Images[{}]
]

2.4.3. Cycles to ...

1) Cycles to Cycles:

TranslatePerm[perm:Cycles[___List],Cycles]:=Cycles@@First[System`Cycles[List@@perm]];

test[TranslatePerm,
TranslatePerm[Cycles[{3,2},{5,6,4}],Cycles],
Cycles[{2,3},{4,5,6}]
]

2) Cycles to Rules:

TranslatePerm[Cycles[cycs___List],Rules]:=Inner[Rule,Flatten[{cycs}],Flatten[RotateLeft/@{cycs}],Rules];

test[TranslatePerm,
TranslatePerm[Cycles[{1,3,2},{4,5}],Rules],
Rules[1->3,3->2,2->1,4->5,5->4]
]

test[TranslatePerm,
TranslatePerm[Cycles[{1,2}],Rules],
Rules[1->2,2->1]
]

Trace[TranslatePerm[Cycles[{102,101}],Rules]]

3) Cycles to Perm:

TranslatePerm[Cycles[cycs___List],{Perm,length_Integer}]:=Perm@PermList[System`Cycles[Reverse/@{cycs}],length]

test[TranslatePerm,
TranslatePerm[Cycles[{1,2,3},{4,7}],{Perm,9}],
Perm[{3,1,2,7,5,6,4,8,9}]
]

test[TranslatePerm,
TranslatePerm[Cycles[{1,2,3},{4,7}],Perm],
Perm[{3,1,2,7,5,6,4}]
]

4) Cycles to Images:

TranslatePerm[perm:Cycles[cycs___List],{Images,length_Integer}]:=Images@PermList[System`Cycles[{cycs}],length];

test[TranslatePerm,
TranslatePerm[Cycles[{1,2,3},{4,7}],Images],
Images[{2,3,1,7,5,6,4}]
]

test[TranslatePerm,
TranslatePerm[-Cycles[{1,2,3},{4,7}],{Images,8}],
-Images[{2,3,1,7,5,6,4,8}]
]

2.4.4. Perm to ...

1) Perm to Perm. Note that if we require a larger length then we add more points:

TranslatePerm[Perm[list_List],{Perm,length_Integer}]:=Perm@PermList[list,length];

2) Perm to Cycles:

TranslatePerm[Perm[list_List],Cycles]:=Cycles@@(Reverse/@First[PermCycles[list]]);

test[TranslatePerm,
TranslatePerm[Perm[{3,4,2,5,1,6,8,7}],Cycles],
Cycles[{5,4,2,3,1},{8,7}]
]

3) Perm to Rules (intermediate translation to Cycles):

TranslatePerm[perm:Perm[_List],Rules]:=TranslatePerm[TranslatePerm[perm,Cycles],Rules];

test[TranslatePerm,
TranslatePerm[Perm[{3,4,2,5,1,6,8,7}],Rules],
Rules[5->4,4->2,2->3,3->1,1->5,8->7,7->8]
]

test[TranslatePerm,
TranslatePerm[Perm[{}],Rules],
Rules[]
]

4) Perm To Images. Again, add points if we require a larger length:

TranslatePerm[Perm[list_List],{Images,length_Integer}]:=Images[Take[Ordering[list]~Join~Range[Length[list]+1,length],length]];

test[TranslatePerm,
TranslatePerm[Perm[{4,2,5,3,1}],Images],
Images[{5,2,4,1,3}]
]

test[TranslatePerm,
TranslatePerm[Perm[{4,2,5,3,1,6,7,8}],Images],
Images[{5,2,4,1,3,6,7,8}]
]

test[TranslatePerm,
TranslatePerm[Perm[{4,2,5,3,1,6,7,8}],{Images,5}],
Images[{5,2,4,1,3}]
]

test[TranslatePerm,
TranslatePerm[Perm[{4,2,5,3,1,6,7,8}],{Images,10}],
Images[{5,2,4,1,3,6,7,8,9,10}]
]

2.4.5. Rules to ...

1) Rules to Rules:

TranslatePerm[perm:Rules[___Rule],Rules]:=perm;

2) Rules to Cycles. Rewriting-rules coding:

m1c[{n_,n_},_]:=Sequence[];
m1c[{n_,other__,n_},_]:={n,other};
m1c[{other___,n_},rules_]:=m1c[{other,n,n/.rules},rules];
mcs[{n_,other___},rules_,{cycs___}]:=With[{cyc=m1c[{n},rules]},mcs[Complement[{other},cyc],rules,{cycs,cyc}]];
mcs[{},_,cycs_]:=cycs;
TranslatePerm[Rules[rules___Rule],Cycles]:=Cycles@@mcs[Range[PermDeg[Rules[rules]]],{rules},{}];

It would be nice to write a one-liner, but I don't know how to do it.

test[TranslatePerm,
TranslatePerm[Rules[1->3,3->1],Cycles],
Cycles[{1,3}]
]

test[TranslatePerm,
TranslatePerm[Rules[1->3,3->2,2->1,7->8,8->7],Cycles],
Cycles[{1,3,2},{7,8}]
]

test[TranslatePerm,
TranslatePerm[-Rules[],Cycles],
-Cycles[]
]

3) Rules to Perm:

TranslatePerm[Rules[rules___Rule],{Perm,length_Integer}]:=Perm@ReplacePart[Range@length,First/@{rules},({#[[2]]}&/@{rules})/.{}->{{}},Partition[Range[Length[{rules}]],1]/.{}->{{}}]

test[TranslatePerm,
TranslatePerm[Rules[1->2,2->5,5->1],{Perm,7}],
Perm[{5,1,3,4,2,6,7}]
]

test[TranslatePerm,
TranslatePerm[-Rules[1->2,2->5,5->1,7->8,8->7],Perm],
-Perm[{5,1,3,4,2,6,8,7}]
]

test[TranslatePerm,
TranslatePerm[Rules[],Perm],
Perm[{}]
]

4) Rules to Images:

TranslatePerm[perm:Rules[___Rule],{Images,length_Integer}]:=Images[Range@length/.List@@perm]

test[TranslatePerm,
TranslatePerm[Rules[1->3,2->4,3->2,4->1],Images],
Images[{3,4,2,1}]
]

2.4.6. Images to ...

1) Images to Images:

TranslatePerm[Images[list_List],{Images,length_Integer}]:=Images[Take[list~Join~Range[Length[list]+1,length],length]];

2) Images to Perm:

TranslatePerm[Images[list_List],{Perm,length_}]:=Perm[Take[Ordering[list]~Join~Range[Length[list]+1,length],length]];

test[TranslatePerm,
TranslatePerm[-Perm[{3,4,2,1,7,6,5}],{Images,9}],
-Images[{4,3,1,2,7,6,5,8,9}]
]

test[TranslatePerm,
TranslatePerm[-Images[{4,3,1,2,7,6,5,8,9}],Perm],
-Perm[{3,4,2,1,7,6,5,8,9}]
]

3) Images to Cycles:

TranslatePerm[Images[list_List],Cycles]:=Cycles@@First[PermCycles[list]];

test[TranslatePerm,
TranslatePerm[-Images[{4,3,1,2,7,6,5}],Cycles],
-Cycles[{1,4,2,3},{5,7}]
]

4) Images to Rules:

TranslatePerm[Images[perm_List],Rules]:=Inner[Rule,Sort[perm],perm,Rules];

test[TranslatePerm,
TranslatePerm[-Images[{4,3,1,2,7,6,5}],Rules],
-Rules[1->4,2->3,3->1,4->2,5->7,7->5]
]

2.4.7. Protection

TranslatePerm[perm_,Cycles|Rules|Perm|Images|{Perm,_Integer}|{Images,_Integer}]:=Throw[Message[TranslatePerm::invalid,perm,"permutation"]];
TranslatePerm[_,notation_]:=Throw[Message[TranslatePerm::unknown,"notation",notation]];
SetNumberOfArguments[TranslatePerm,2];
Protect[TranslatePerm];

Catch@TranslatePerm[Hello,Cycles]

Catch@TranslatePerm[Hello,GoodBye]

And@@alltests[TranslatePerm]

2.5. Point images. OnPoints

2.5.1. Definitions and examples

This is usually treated as a Power function on permutations. I don't like that at all.

?OnPoints

OnPoints is listable on its first argument:

OnPoints[list_List,perm_]:=OnPoints[#,perm]&/@list;

Permutations acting on points:

OnPoints[p_Integer,ID]:=p;

HACK: OnPoints with negative permutations: forget about sign. Onpoints is not defined on sums or multiplications for scalars other than -1.

(* Sign points *)
OnPoints[0,-perm_]:=-1;
OnPoints[-1,-perm_]:=0;
OnPoints[0,perm_]:=0;
OnPoints[-1,perm_]:=-1;
(* Other points *)
OnPoints[p_Integer,-perm_]:=OnPoints[p,perm];

1) Perm notation:

OnPoints[p_Integer,Perm[list:{___,p_,___}]]:=Position[list,p][[1,1]];
OnPoints[p_Integer,Perm[{___}]]:=p;

test[OnPoints,
OnPoints[5,Perm[{3,4,2,1}]],
5
]

test[OnPoints,
OnPoints[4,Perm[{3,4,2,1}]],
2
]

test[OnPoints,
OnPoints[1,Perm[{}]],
1
]

The function OnPoints is reasonably fast in Perm notation

Trace[OnPoints[4,Perm[{3,5,4,2,1}]]]

2) Rules notation:

OnPoints[p_Integer,Rules[rules___Rule]]:=p/.{rules};

test[OnPoints,
OnPoints[3,Rules[1->3,3->5,5->1,6->7,7->6]],
5
]

test[OnPoints,
OnPoints[3,-Rules[]],
3
]

test[OnPoints,
OnPoints[-1,-Rules[]],
0
]

Trace[OnPoints[3,Rules[1->3,3->5,5->1,6->7,7->6]]]

3) Cycles notation:

OnPoints[p_Integer,Cycles[___List,{___,p_,q_,___},___List]]:=q;
OnPoints[p_Integer,Cycles[___List,{q_,___,p_},___List]]:=q;
OnPoints[p_Integer,Cycles[___List]]:=p;

test[OnPoints,
OnPoints[3,Cycles[{1,5},{2,3,4}]],
4
]

test[OnPoints,
OnPoints[3,Cycles[{1,5},{2,4,3}]],
2
]

test[OnPoints,
OnPoints[3,Cycles[{1,2}]],
3]

test[OnPoints,
OnPoints[3,Cycles[]],
3
]

test[OnPoints,
OnPoints[{1,2,3,4,5},Cycles[{1,5,2,4,3}]],
{5,4,1,3,2}
]

test[OnPoints,
OnPoints[1,-Cycles[{2,3}]],
1
]

test[OnPoints,
OnPoints[2,-Cycles[{2,3}]],
3
]

4) Images notation:

OnPoints[p_Integer,Images[list:{___,p_,___}]]:=list[[p]];
OnPoints[p_Integer,Images[{___}]]:=p;

test[OnPoints,
OnPoints[1,Images[{2,3,4,1}]],
2
]

test[OnPoints,
OnPoints[5,Images[{2,3,4,1}]],
5
]

Protection:

SetNumberOfArguments[OnPoints,2];
Protect[OnPoints];

It is not defined on linear combinations:

OnPoints[1,Perm[{1,2}]+Perm[{2,1}]]

OnPoints[1,3Perm[{1,2}]]

2.5.2. Timings

First we compare different notations:

Timing[p=1;Do[p=OnPoints[p,Perm[{5,1,6,2,3,4}]],{i,1,10000}];p]

Timing[p=1;Do[p=OnPoints[p,Images[{2,4,5,6,1,3}]],{i,1,10000}];p]

Timing[p=1;Do[p=OnPoints[p,Rules[1->2,2->4,4->6,6->3,3->5,5->1]],{i,1,10000}];p]

Timing[p=1;Do[p=OnPoints[p,Cycles[{5,1,2,4,6,3}]],{i,1,10000}];p]

Now we include irrelevant points:

Timing[p=1;Do[p=OnPoints[p,Perm[{5,1,6,2,3,4,7,8,9,10,11,12}]],{i,1,10000}];p]

Timing[p=1;Do[p=OnPoints[p,Images[{2,4,5,6,1,3,7,8,9,10,11,12}]],{i,1,10000}];p]

Timing[p=1;Do[p=OnPoints[p,Rules[1->2,2->4,4->6,6->3,3->5,5->1,7->7,8->8,9->9,10->10,11->11,12->12]],{i,1,10000}];p]

Timing[p=1;Do[p=OnPoints[p,Cycles[{5,1,2,4,6,3},{7},{8},{9},{10},{11},{12}]],{i,1,10000}];p]

It is clear that a lot of time is lost eliminating the points at the end of the permutations.

2.6. Composition of permutations. PermProduct

Composition of two permutations is defined. Composition of permutations in different notations is not allowed. Note that this is an important point to speed up calculations. In particular, this is why we want to ensure that all permutations in Perm or Images notation have the same length.

Linear combinations (we shall need an algebra to work with Young tableaux):

PermProduct[perm1_,expr_Plus]:=PermProduct[perm1,#]&/@expr;
PermProduct[expr_Plus,perm2_]:=PermProduct[#,perm2]&/@expr;
PermProduct[x_?NumericQ perm1_,perm2_]:=x PermProduct[perm1,perm2];
PermProduct[perm1_,y_?NumericQ perm2_]:=y PermProduct[perm1,perm2];
PermProduct[perm_,0]:=0;
PermProduct[0,perm_]:=0;

Identity:

PermProduct[]:=ID;
PermProduct[perm_,ID]:=perm;
PermProduct[ID,perm_]:=perm;

OneIdentity (this is given to be consistent with PermPower[perm, 1] below):

PermProduct[perm_]:=perm;

Associativity (it doesn't matter which side we choose):

PermProduct[perm1_,perm2_,perm3__]:=PermProduct[perm1,PermProduct[perm2,perm3]];

Thread over lists:

PermProduct[list_List,perm_]:=PermProduct[#,perm]&/@list;
PermProduct[perm_,list_List]:=PermProduct[perm,#]&/@list;

Linearity:

PermProduct[Perm[{1,2,3}]+3Perm[{2,1,3}],Perm[{3,2,1}]]

1) Perm notation (note that we do not check lengths). This definition essentially defines the convention of multiplications: from left to right:

PermProduct[Perm[list_List],Perm[{}]]:=Perm[list];
PermProduct[Perm[{}],Perm[list_List]]:=Perm[list];
PermProduct[Perm[list1_List],Perm[list2_List]]:=Perm[list1[[list2]]];

test[PermProduct,
PermProduct[Perm[{3,4,5,1,2,6}],Perm[{5,3,4,1,6,2}]],
Perm[{2,5,1,3,6,4}]
]

This is rather efficient:

Trace[PermProduct[Perm[{3,4,5,1,2,6}],Perm[{5,3,4,1,6,2}]]]

2) Rules notation:

PermProduct[Rules[r1___Rule],Rules[r2___Rule]]:=Inner[Rule,#,#/.{r1}/.{r2},Rules]&@Union[Level[{r1,r2},{-1}]]

test[PermProduct,
PermProduct[Rules[1->4,4->2,2->5,5->3,3->1],Rules[1->4,4->3,3->2,2->6,6->5,5->1]],
Rules[1->3,2->1,3->4,4->6,5->2,6->5]
]

test[PermProduct,
PermProduct[-Rules[1->4,4->2,2->5,5->3,3->1],-Rules[1->4,4->3,3->2,2->6,6->5,5->1]],
Rules[1->3,2->1,3->4,4->6,5->2,6->5]
]

test[PermProduct,
PermProduct[Rules[],-Rules[]],
-Rules[]
]

This is not very efficient:

Trace[PermProduct[Rules[1->4,4->2,2->5,5->3,3->1],Rules[1->4,4->3,3->2,2->6,6->5,5->1]]]

3) Cycles notation (just translation to Rules). This should be improved:

PermProduct[Cycles[],perm:Cycles[___List]]:=perm;
PermProduct[perm:Cycles[___List],Cycles[]]:=perm;
PermProduct[Cycles[cycs1___],Cycles[cycs2___]]:=Apply[Cycles,First@PermutationProduct[System`Cycles[{cycs1}],System`Cycles[{cycs2}]]];

test[PermProduct,
PermProduct[Cycles[{1,3}],Cycles[{1,2},{3,4}]],
Cycles[{1,4,3,2}]
]

test[PermProduct,
PermProduct[Cycles[{1,5,3}],Cycles[{1,2,3,4}]],
Cycles[{1,5,4},{2,3}]
]

4) Images notation:

PermProduct[Images[list1_List],Images[list2_List]]:=Images[Union[list1,list2]/.Inner[Rule,Sort[list1],list1,List]/.Inner[Rule,Sort[list2],list2,List]];

This definition is insensitive to the lengths of the lists of points of the permutations.

test[PermProduct,
PermProduct[Images[{2,3,4,1,5,6}],-Images[{3,4,5,1,6,2}]],
-Images[{4,5,1,3,6,2}]
]

Trace[PermProduct[Images[{2,3,4,1,5,6}],Images[{3,4,5,1,6,2}]]]

This is a very important point, because there is a choice between two possibilities. Our choice means that given a list omega we perform permutations in the form omega ∘ perm = omega[[ perm ]]. For example, in Perm notation, if we are given two permutations perm1 and perm2, we can multiply them as perm1 ∘ perm2, which means omega ∘ ( perm1 ∘ perm2 ) = omega [[ perm1 ∘ perm2 ]] = omega[[ perm1 ]][[ perm2 ]] = omega[[ perm1[[ perm2 ]] ]]. Therefore it is consistent to define perm1 ∘ perm2 = perm1[[ perm2 ]]. This is good to remember the order of everything. A posteriori, the main motivation is having, given a point alpha, that      ( alpha^perm1 )^perm2 = alpha^( perm1.perm2)

test[PermProduct,
PermProduct[Perm[{2,3,1,4}],Perm[{2,1,4,3}]],
Perm[{3,2,4,1}]
]

test[PermProduct,
PermProduct[Perm[{2,1,4,3}],Perm[{2,3,1,4}]],
Perm[{1,4,2,3}]
]

This is the choice in GAP, the choice in Butler's book, and the choice in Renato's papers.

This is the example given in GAP Reference, chapter 40:

test[PermProduct,
PermProduct[Cycles[{1,2,3}],Cycles[{2,3,4}]],
Cycles[{1,3},{2,4}]
]

test[OnPoints,
OnPoints[17,Cycles[{2,5,17,9,8}]],
9
]

test[PermProduct,
PermProduct[Cycles[{2,3,4}],Cycles[{1,2,3}]],
Cycles[{1,2},{3,4}]
]

Note that in the Perm and Images notations we do (implicitly) check the lengths, in the sense that errors quickly crop up:

test[PermProduct,
PermProduct[Perm[{2,4,1,3}],Perm[{2,1,3}]],
Perm[{4,2,1}]
]

test[PermProduct,
PermProduct[Perm[{2,1,3}],Perm[{2,4,1,3}]],
$Failed,
True]

Cosets. We first check if the permutation is already in the set or not:

PermProduct[Group[perms___],perm_]:=If[Or@@(PermEqual[#,perm]&/@{perms}),Group[perms],PermProduct[#,perm]&/@Coset[perms]];
PermProduct[perm_,Group[perms___]]:=If[Or@@(PermEqual[perm,#]&/@{perms}),Group[perms],PermProduct[perm,#]&/@Coset[perms]];

Protection. PermProduct can have any number of arguments:

Protection[PermProduct];

2.7. Inverse permutation. InversePerm

2.7.1. Definitions

The only nontrivial part is the use of the function Ordering:

?InversePerm

InversePerm[Perm[list_List]]:=Perm[InversePermutation[list]];
InversePerm[Images[list_List]]:=Images[InversePermutation[list]];
InversePerm[Cycles[cycs___List]]:=Apply[Cycles,First@InversePermutation[System`Cycles[{cycs}]]];
InversePerm[Rules[rules___Rule]]:=Reverse/@Rules[rules];
InversePerm[ID]:=ID;
InversePerm[x_?NumericQ perm_]:=1/x InversePerm[perm];
SetNumberOfArguments[InversePerm,1];
Protect[InversePerm];

2.7.2. Examples

test[InversePerm,
InversePerm[Perm[{2,3,4,1}]],
Perm[{4,1,2,3}]
]

test[InversePerm,
InversePerm[Perm[{}]],
Perm[{}]
]

test[InversePerm,
InversePerm[Images[{4,1,2,3}]],
Images[{2,3,4,1}]
]

test[InversePerm,
InversePerm[-Images[{1,2}]],
-Images[{1,2}]
]

test[InversePerm,
InversePerm[-Cycles[{2,3},{5,6,1}]],
-Cycles[{1,6,5},{2,3}]
]

test[InversePerm,
InversePerm[-Cycles[]],
-Cycles[]
]

test[InversePerm,
InversePerm[Rules[1->3,3->5,5->1]],
Rules[3->1,5->3,1->5]
]

test[InversePerm,
InversePerm[Rules[]],
Rules[]
]

InversePerm is not very efficient:

Trace[InversePerm[Perm[{2,3,4,1}]]]

Trace[InversePerm[Rules[1->3,3->5,5->1]]]

2.8. Permute lists. PermuteList, PermutationFromTo

This function permutes the list according to permutation perm. We neglect signs. This function is used twice in xTensor and then here we use the (private) function SignedPermuteList (see section 3.6 below).

?PermuteList

PermuteList[list_,-perm_]:=PermuteList[list,perm];
PermuteList[list_,perm_]:=list[[First@TranslatePerm[perm,{Perm,Length[list]}]]];
SetNumberOfArguments[PermuteList,2];
Protect[PermuteList];

test[PermuteList,
PermuteList[{1,3,4,2},Images[{3,4,2,1}]],
{2,4,1,3}
]

test[PermuteList,
PermuteList[head[2,4,1,3],Images[{3,4,2,1}]],
head[3,1,2,4]
]

PermuteList[head[],ID]//Trace

We introduce a definition that always works, even if the elements are repeated, in the sense that it always returns a valid permutation

?PermutationFromTo

PermutationFromTo::diff="Lists`1` and `2` contain different points";
PermutationFromTo[list1_List,list2_List]:=PermProduct[InversePerm@Images@Ordering[list1],Images@Ordering[list2]]/;Sort[list1]===Sort[list2];
PermutationFromTo[list1:head_[___],list2:head_[___]]:=PermutationFromTo[List@@list1,List@@list2];
PermutationFromTo[list1_List,list2_List]:=Throw[Message[PermutationFromTo::diff,list1,list2]]/;Sort[list1]=!=Sort[list2];
SetNumberOfArguments[PermutationFromTo,2];
Protect[PermutationFromTo];

test[PermutationFromTo,
PermutationFromTo[{a,b,-c,b,a,d},{b,a,-c,d,a,b}],
Images[{2,1,3,6,5,4}]
]

test[PermuteList,
PermuteList[{a,b,-c,b,a,d},Images[{2,1,3,6,5,4}]],
{b,a,-c,d,a,b}
]

The internal order of the repeated points is not changed:

test[PermuteList,
PermuteList[{a1,b1,-c,b2,a2,d},Images[{2,1,3,6,5,4}]],
{b1,a1,-c,d,a2,b2}
]

test[PermutationFromTo,
PermutationFromTo[{1,3,4,2},{2,4,1,3}],
Images[{3,4,2,1}]
]

test[PermutationFromTo,
PermutationFromTo[{1,3,4,2},{2,4,1,3}],
Images[{3,4,2,1}]
]

Catch@PermutationFromTo[{1,3,4,2},{2,4,1,2}]

2.9. RandomPerm

RandomPerm gives a random permutation of given degree and notation. It is constructed in Perm notation and then translated. Note that the given permutation could have lower degree than that asked for. It would be highly counterintuitive to impose that the last element must be necessarily moved, but not the others.

?RandomPerm

rand[n_]:=Random[Integer,{1,n}];
addrand[n_][list_List,m_]:=Append[list,Complement[Range[n],list][[rand[m]]]];
RandomPerm[deg_]:=RandomPerm[deg,Perm];
RandomPerm[deg_,Perm]:=Perm[Fold[addrand[deg],{},Reverse@Range[deg]]];
RandomPerm[deg_,notation_]:=TranslatePerm[RandomPerm[deg,Perm],notation];
SetNumberOfArguments[RandomPerm,{1,2}];
Protect[RandomPerm];

RandomPerm[15,Rules]

RandomPerm[4,{Images,7}]

With notations Rules and Cycles and degree 1 we can only get the identity:

RandomPerm[0,Rules]

RandomPerm[1,Cycles]

2.10. PermSignature

Signature of permutation. We consider the sign as a transposition more.

PermSignature[Perm[list_List]]:=Signature[list];
PermSignature[Images[list_List]]:=Signature[list];
signofcycle[cyc_List]:=-(-1)^Length[cyc];
PermSignature[perm_Cycles]:=Apply[Times,signofcycle/@perm];
PermSignature[perm_Rules]:=PermSignature[TranslatePerm[perm,Cycles]];
PermSignature[-perm_]:=-PermSignature[perm];
SetNumberOfArguments[PermSignature,1];
Protect[PermSignature];

3. Permutation Defs. Notation independent

(************* 3. Permutations. Notation independent *************)

{TimeUsed[],SessionTime[],TimeUsed[]/SessionTime[]}

3.0. Comments

This section defines all those actions on permutations whose structure is the same on different notations for the permutations.

If[$ReadingVerbose,Print["Reading section 3: Permutations. Notation independent"],Null,Null]

3.1. Sorting with a base

Given a set of points Ω, we can define inequivalent ways to sort them. This is very important because permutations will be sorted according to the order of points. Different orderings will be given using "bases": ordered subsets of Ω.

Note that we do not introduce a head Base because List is natural enough.

We define three different functions to work with orderings. They have independent codings, but they could be written in terms of just one of them:

1) First we define the function that sorts elements according to a given ordered set B. Note that points not appearing in the B list are sorted last, and sorted according to Sort:

?SortB

SortB[list_List,B_List]:=DeleteCases[B,Alternatives@@Complement[B,list]]~Join~Complement[list,B];
SortB[list:{_},_]:=list;
SortB[{},_]:={};
SortB[list_List]:=Sort[list];
SetNumberOfArguments[SortB,2];
Protect[SortB];

test[SortB,
SortB[{6,1,4,5,2},{4,2,3,5}],
{4,2,5,1,6}
]

SortB[{6,1,4,5,2},{4,2,3,5}]//Trace

2) This rather efficient function finds the minimum of a list with respect to a base:

?MinB

MinB[{___,p_,___},{p_,___}]:=p;
MinB[list_List,{_,x___}]:=MinB[list,{x}];
MinB[list_List,{}]:=Min[list];
MinB[list_List]:=Min[list];
SetNumberOfArguments[MinB,2];
Protect[MinB];

MinB[{6,1,4,5,2},{4,2,3,5}]//Trace

test[MinB,
MinB[{6,1,4,5,2},{4,2,3,5}],
4
]

MinB[{6,1,4,5,2},{3,7}]//Trace

test[MinB,
MinB[{6,1,4,5,2},{3,7}],
1
]

test[SortB,
SortB[{6,1,4,5,2},{3,7}],
{1,2,4,5,6}
]

3) We assume that the base does not have repeated points. If one of the points is in the base but not the other, then the canonical order is having the former point first. If none of the points is in the base then we use normal lexicographical order. This function is private.

LessB[{x_,y_,z__},base_List]:=LessB[{x,y},base]&&LessB[{y,z},base];
LessB[{x_,x_},_List]:=False;
LessB[{x_,y_},{___,x_,___,y_,___}]:=True;
LessB[{x_,y_},{___,y_,___,x_,___}]:=False;
LessB[{x_,_},{___,x_,___}]:=True;
LessB[{_,x_},{___,x_,___}]:=False;
LessB[{x_,y_},_List]:=Less[x,y];

test[LessB,
LessB[{4,2,1},{4,2,3,5}],
True
]

test[LessB,
LessB[{1,2,4},{4,2,3,5}],
False
]

Less[{1,2,4}]

test[SortB,
SortB[{1,2,4},{4,2,3,5}],
{4,2,1}
]

3.2. Sorting permutations. PermOrderedQ, PermSort

3.2.0. Comments

There are two main ways of sorting permutations: using images of points or using the Perm notation. It turns out that using the former for a set of permutations is equivalent to using the latter for the set of inverses. PermOrderedQ returns Null when both permutations are equal (that is, they give same images for same points; this differs from the usual behaviour of OrderedQ).

3.2.1. PermOrderedQ and PermSort

Note that a linear combination of permutations would give errors because OnPoints is not defined on linear combinations. In PermOrderedQ head is needed but irrelevant. This definition means that we are sorting according to Images, and not Perm.

?PermOrderedQ

PermOrderedQ[head_[perm1_/;Head[perm1]=!=Times,-perm2_]]:=True;
PermOrderedQ[head_[-perm1_,perm2_/;Head[perm2]=!=Times]]:=False;
PermOrderedQ[head_[-perm1_,-perm2_]]:=PermOrderedQ[head[perm1,perm2]];

PermOrderedQ[head_[perm1_,perm2_]]:=Module[{len=Max[PermLength@perm1,PermLength@perm2],i,p1,p2},
For[i=1,i≤len,++i,
p1=OnPoints[i,perm1];
p2=OnPoints[i,perm2];
Which[
p1<p2,Return[True],
p1>p2,Return[False]]
];
Null
];

We define a new function that depends on a base of points.

This is the definition given by Renato for two different permutations. It is not equivalent to my definition below:

PermOrderedQ[head_[perm1_,perm2_],base_List]:=Module[{deg=Max[PermDeg@perm1,PermDeg@perm2],i,p1,p2},
For[i=1,i≤Length[base],++i,
p1=OnPoints[base[[i]],perm1];If[FreeQ[base,p1],p1=p1+deg];
p2=OnPoints[base[[i]],perm2];If[FreeQ[base,p2],p2=p2+deg];
Which[
p1<p2,Return[True],
p1>p2,Return[False]]
];
Null
];

(The addition of n is equivalent to our completion of the base.) The essential difference is that we use LessB comparisons and not the usual Less comparison.

PermOrderedQ[head_[perm1_,perm2_],base_List]:=Module[{len=Max[PermLength@perm1,PermLength@perm2],i,b,p1,p2,base2},
base2=Join[base,Complement[Range[len],base]];
For[i=1,i≤len,++i,b=base2[[i]];
p1=OnPoints[b,perm1];
p2=OnPoints[b,perm2];
Which[
LessB[{p1,p2},base],Return[True],
LessB[{p2,p1},base],Return[False]]
];
Null
];

Corresponding Sort function:

PermSort[head_[perms___]]:=head@@Sort[{perms},PermOrderedQ[{#1,#2}]&];
PermSort[head_[perms___],base_List]:=head@@Sort[{perms},PermOrderedQ[{#1,#2},base]&];

Protection:

SetNumberOfArguments[PermOrderedQ,{1,2}];
SetNumberOfArguments[PermSort,{1,2}];
Protect[PermOrderedQ,PermSort];

These permutations are not sorted lexicographically:

PermSort[{Perm[{1,2,3}],Perm[{1,3,2}],Perm[{2,1,3}],Perm[{2,3,1}],Perm[{3,1,2}],Perm[{3,2,1}]}]

but are sorted in such a way that their inverses are in lexicographical order:

InversePerm/@%

The problem is in these two permutations:

test[OnPoints,
OnPoints[{1,2,3},Perm[{3,1,2}]],
{2,3,1}
]

test[OnPoints,
OnPoints[{1,2,3},Perm[{2,3,1}]],
{3,1,2}
]

Sorting the identity. Different forms of the identity are not sorted in a unique way. Compare the relative positions of Cycles[] and Perm[{}] in these two lists:

PermSort[{Perm[{2,1,3}],ID,Cycles[],Perm[{}],ID}]

PermSort[{Perm[{2,1,3}],ID,Perm[{}],Cycles[],Perm[{}],ID}]

3.2.2. PermEqual, PermLess, PermGreater, PermLessEqual, PermGreaterEqual

We generalize comparatives (not for B-sorting):

PermEqual[perms__]:=And@@(#===Null&/@PermOrderedQ/@Partition[{perms},2,1]);
PermLessEqual[perms__]:={perms}===PermSort[{perms}];
PermGreaterEqual[perms__]:=PermLessEqual@@Reverse[{perms}];
PermLess[perms__]:=And@@(PermOrderedQ/@Partition[{perms},2,1]/.Null->False);
PermGreater[perms__]:=PermLess@@Reverse[{perms}];
Protect[PermEqual,PermLess,PermGreater,PermLessEqual,PermGreaterEqual];

Examples:

test[PermEqual,
PermEqual[Cycles[],Perm[{1,2,3,4}]],
True
]

test[PermEqual,
PermEqual[Cycles[],ID],
True
]

test[PermEqual,
PermEqual[Perm[{}],Perm[{1,2,3,4}]],
True
]

test[PermEqual,
PermEqual[Perm[{}],Perm[{1,2,3,4}]],
True
]

test[PermEqual,
PermLess[Cycles[],Perm[{1,2,3,4}]],
False
]

test[PermEqual,
PermLessEqual[ID,Perm[{1,2,3,4}]],
True
]

test[PermEqual,
PermLess[Perm[{2,3,1}],Perm[{3,1,2}]],
False
]

Examples in GAP Reference, chapter 40 (the first should give True and the second False):

test[PermEqual,
PermLess[Cycles[{1,2,3}],Cycles[{1,3,2}]],
True
]

test[PermEqual,
PermLess[Cycles[{1,3,2,4}],Cycles[{1,3,4,2}]],
False
]

3.3. Set of generators. GenSet

Generator sets will be given as lists of permutations, with flat head GenSet. By convention we assume that we will never have the identity on a generator set, but we do not check that.

?GenSet

GenSet[x___,GenSet[y___],z___]:=GenSet[x,y,z];

GenSet extracts the generating set from a strong generating set:

GenSet[StrongGenSet[_,GS_GenSet]]:=GS;

Special cases:

GenSet[(Symmetric|Antisymmetric)[{}|{_}]]:=GenSet[];
GenSet[Symmetric[list_List]]:=DeleteDuplicates@GenSet[Cycles[list[[{1,2}]]],Cycles[list]];
GenSet[Antisymmetric[list_List]]:=DeleteDuplicates@GenSet[-Cycles[list[[{1,2}]]],-(-1)^Length[list]Cycles[list]];

Empty generating sets represents the group Group[ ID ] and therefore it is reasonable to assign them degree 0:

PermDeg[GenSet[]]

Protect[GenSet];

3.4. Powers. PermPower

?PermPower

For powers of index larger than 10 we switch to a recursive algorithm, far more efficient than a simple listing of permutations:

PermPower[perm_,0]:=ID[perm];
PermPower[perm_,1]:=perm;
PermPower[perm_,n_Integer?Negative]:=PermPower[InversePerm[perm],-n];
PermPower[perm_,n_Integer?Positive]:=PermProduct@@Table[perm,{n}]/;n<10;
PermPower[perm_,n_Integer?Positive]:=With[{m=IntegerPart[n/2]},PermProduct[PermPower[perm,m],PermPower[perm,n-m]]];
PermPower[perm_,n_Integer?Positive]:=Module[{pp,p},
pp[p_,0]:=IP[p];
pp[p_,1]:=p;
pp[p_,k_]:=pp[p,k]=With[{m=IntegerPart[k/2]},PermProduct[pp[p,m],pp[p,k-m]]];
pp[perm,n]
];
SetNumberOfArguments[PermPower,2];
Protect[PermPower];

Examples:

test[PermPower,
PermPower[Perm[{}],7],
Perm[{}]
]

test[PermPower,
PermPower[ID,7],
ID
]

perm=Perm[{5,4,3,2,1}];

test[PermPower,
PermPower[perm,0],
Perm[{1,2,3,4,5}]
]

This is a permutation of order 4:

AbsoluteTiming[PermPower[perm,4]]

AbsoluteTiming[PermPower[perm,40]]

AbsoluteTiming[PermPower[perm,4000]]

AbsoluteTiming[PermPower[perm,4000000]]

AbsoluteTiming[PermPower[perm,4567891]]

AbsoluteTiming[PermPower[perm,-4567891]]

The growth is clearly not linear!

Note that for negative n InversePerm is used  just once.

3.5. Dimino algorithm. Dimino, Group

And finally, this the simple Dimino algorithm to generate a group by generators. MakeGroup enlarges subgroup (assumed to be generated by all but the last generators in GS; not checked) to the group generated by all generators in GS. Note that we use MemberQ and FreeQ to test for membership of permutations; that is only valid for {Perm,length} and {Images,length} notations and hence we have to change to Perm notation initially.

This is algorithm 3.4 of Butler.

cyclicGroup[perm_]:=With[{id=ID[perm]},Group@@RotateRight[NestWhileList[PermProduct[#,perm]&,perm,(#=!=id)&]]];
RightCoset[subgroup_Group,perm_]:=Map[PermProduct[#,perm]&,subgroup];
MakeGroup[H_Group,gs_GenSet]:=If[MemberQ[H,Last[gs]],H,Module[{
ordH=Length[H],
elts=H~Join~RightCoset[H,Last[gs]],
ncoset=2},
While[ordH ncoset≤Length[elts],
Map[If[FreeQ[elts,#],elts=elts~Join~RightCoset[H,#]]&,Map[PermProduct[elts[[ordH (ncoset-1)+1]],#]&,gs]];
++ncoset
];
elts]];

?Dimino

Dimino[StrongGenSet[base_,GS_GenSet]]:=Dimino[GS];
Dimino[GS_GenSet]:=Module[{GSPerm=TranslatePerm[GS,{Perm,PermDeg[GS]}]},
TranslatePerm[
Fold[
MakeGroup[#1,#2]&,
cyclicGroup[First[GSPerm]],
Rest[FoldList[Append[#1,#2]&,GenSet[First[GSPerm]],List@@Rest[GSPerm]]]
],NotationOfPerm[GS]]
];
Dimino[GenSet[]]:=Group[ID];

Dimino[Symmetric[list_]]:=Group@@(Images/@Permutations[list]);
Dimino[Antisymmetric[list_]]:=Group@@((Signature[#]Images[#])&/@Permutations[list]);

SetNumberOfArguments[Dimino,1];
Protect[Dimino,Group];

Examples:

test[Dimino,
Dimino[GenSet[Perm[{4,1,2,3}]]],
Group[Perm[{1,2,3,4}],Perm[{4,1,2,3}],Perm[{3,4,1,2}],Perm[{2,3,4,1}]]
]

test[Dimino,
Dimino[GenSet[-Perm[{4,1,2,3}]]],
Group[Perm[{1,2,3,4}],-Perm[{4,1,2,3}],Perm[{3,4,1,2}],-Perm[{2,3,4,1}]]
]

test[Dimino,
Dimino[GenSet[Images[{2,3,4,1}]]],
Group[Images[{1,2,3,4}],Images[{2,3,4,1}],Images[{3,4,1,2}],Images[{4,1,2,3}]]
]

test[Dimino,
Dimino[GenSet[Rules[1->2,2->3,3->4,4->1]]],
Group[Rules[],Rules[2->3,3->4,4->1,1->2],Rules[3->1,1->3,4->2,2->4],Rules[4->3,3->2,2->1,1->4]]]

Dimino[GenSet[Cycles[{1,2,3,4}]]]

Permutations in Cycles notation are not being returned sorted. Can we avoid that?

Timing[Sort[Dimino[GenSet[Perm[{4,1,2,3}],Perm[{2,1,3,4}]]]]===Group@@Sort[Perm/@Permutations[{1,2,3,4}]]]

Recall that

{5!,6!,7!}

Timing[Length@Dimino[GenSet[Perm[{5,1,2,3,4}],Perm[{2,1,3,4,5}]]]]

Timing[Length@Dimino[GenSet[Perm[{6,1,2,3,4,5}],Perm[{2,1,3,4,5,6}]]]]

Timing[Length@Dimino[GenSet[Perm[{7,1,2,3,4,5,6}],Perm[{2,1}]]]]

This is quite slow now... Compare with

Timing[Length@Permutations[Range[7]]]

3.6. Stabilization I

3.6.1. StablePoints

This function gives the points between 1 and n which remain at place under the permutation p, or under a generating set of permutations

?StablePoints

StablePoints[perm_]:=StablePoints[perm,PermLength[perm]];
StablePoints[perm_?PermQ,n_Integer]:=Cases[Transpose[{Range[n],OnPoints[Range[n],perm]}],{x_,x_}->x];
StablePoints[GS:GenSet[__],n_Integer]:=Intersection@@(StablePoints[#,n]&/@GS);
StablePoints[GenSet[],n_Integer]:=Range[n];
SetNumberOfArguments[StablePoints,{1,2}];
Protect[StablePoints];

test[StablePoints,
StablePoints[Perm[{3,2,1,4,7,5,6}]],
{2,4}
]

test[StablePoints,
StablePoints[Perm[{3,2,1,4,7,5,6}],8],
{2,4,8}
]

test[StablePoints,
StablePoints[Cycles[{3,4}],4],
{1,2}
]

test[StablePoints,
StablePoints[ID,0],
{}
]

test[StablePoints,
StablePoints[ID,8],
{1,2,3,4,5,6,7,8}
]

test[StablePoints,
StablePoints[GenSet[Perm[{1,2,3,5,4}],Perm[{1,2,4,5,3}]],6],
{1,2,6}
]

test[StablePoints,
StablePoints[GenSet[Perm[{1,2,3,5,4}],Perm[{1,2,4,5,3}]],5],
{1,2}
]

test[StablePoints,
StablePoints[GenSet[],5],
{1,2,3,4,5}
]

test[StablePoints,
StablePoints[GenSet[],0],
{}
]

And@@alltests[StablePoints]

3.6.2. NonStablePoints

This function gives the minimal number of points (trying to minimize them as well), including initB, such that none of the generators GS fixes all of them.

?NonStablePoints

NonStablePoints[initB_List,GenSet[]]:=initB
NonStablePoints[initB_List,GS:GenSet[__]]:=Module[{
SubsetQ,stablepoints,base=initB,candidates,check,len=PermLength[GS]},
SubsetQ[large_List,small_List]:=And@@(MemberQ[large,#]&)/@small;
stablepoints=List@@(StablePoints[#,len]&/@GS);
candidates=DeleteCases[Join[Complement[Range[len],#],#]&@Flatten[Union/@Sort@Split@Sort@Flatten@stablepoints],Alternatives@@Intersection@@stablepoints];
While[True,
candidates=DeleteCases[candidates,Alternatives@@base];
check=SubsetQ[#,base]&/@stablepoints;
Which[
(* Correct base. Break while loop *)
Not[Or@@check],Break[],
(* No candidate elements. Throw error *)
candidates==={},Throw["Unable to find enough nonstable points."],
(* Incorrect base. Add new point *)
True,For [i=1,i≤Length[check],++i,
If[check[[i]],AppendTo[base,First[Cases[candidates,x_/;FreeQ[stablepoints[[i]],x]]]];Break[]]
]]];
base
];
SetNumberOfArguments[NonStablePoints,2];
Protect[NonStablePoints];

NonStablePoints[{},GenSet[]]

From example 4 below:

a=Cycles[{1,8,9},{2,11,15},{3,10,12},{4,14,19},{5,16,17},{6,21,20},{7,13,18}];
b=Cycles[{9,18,20},{12,19,17}];
c=Cycles[{10,21,11},{13,16,14}];

test[NonStablePoints,
NonStablePoints[{},GenSet[a,b,c]],
{9,10}
]

test[NonStablePoints,
NonStablePoints[{},GenSet[b,c,a]],
{9,10}
]

test[NonStablePoints,
NonStablePoints[{3},GenSet[a,b,c]],
{3,9,10}
]

test[NonStablePoints,
NonStablePoints[{9},GenSet[a,b,c]],
{9,10}
]

3.6.3. Pointwise stabilization: Stabilizer

This function gives the subset of permutations that leave stable the set of points. Note that when there is no stabilizing permutation we return GenSet[], and not the identity permutation, if it was not in GS.

?Stabilizer

Stabilizer[{},GS_GenSet]:=GS;
Stabilizer[{p_Integer},GS_GenSet]:=GenSet@@Cases[{p,OnPoints[p,#],#}&/@GS,{x_,x_,y_}->y];
Stabilizer[{p1_Integer,ps__Integer},GS_GenSet]:=Stabilizer[{ps},Stabilizer[{p1},GS]];
Stabilizer[ps_List,Symmetric[list_]]:=Symmetric[Complement[list,ps]];
Stabilizer[ps_List,Antisymmetric[list_]]:=Antisymmetric[Complement[list,ps]];
SetNumberOfArguments[Stabilizer,2];

The function is not protected because it will be overloaded below with definitions on SGSs.

Stabilizer[{1},GenSet[Perm[{3,2,1,4,7,5,6}],Perm[{1,2,4,3,6,5,7}]]]

Stabilizer[{1,3},GenSet[Perm[{3,2,1,4,7,5,6}],Perm[{1,2,4,3,6,5,7}]]]

Remember that an empty generating set actually represents the Group[ ID ].

There is an important point to consider: The stabilizer of the generating set of a group is not necessarily a generating set of the stabilizer subgroup. For example

Stabilizer[{2},GenSet[Perm[{2,1,3}],Perm[{1,3,2}]]]

and not GenSet[Perm[{3,2,1}]]. This problem will be solved below using strong generating sets.

3.6.4. Set stabilization: SetStabilizer

Now we do not require that all points in a given list are stable under the stabilizer. We just require that those points do not leave the list.

?SetStabilizer

SetStabilizer[{},GS_GenSet]:=GS;
SetStabilizer[{p_Integer},GS_GenSet]:=Stabilizer[{p},GS];
SetStabilizer[list_List,GS_GenSet]:=setstabilizer[list,GS,GenSet[]];
setstabilizer[list_List,GenSet[perm1_,perms___],GenSet[stperms___]]:=setstabilizer[list,GenSet[perms],If[Complement[OnPoints[list,perm1],list]==={},GenSet[stperms,perm1],GenSet[stperms]]];
setstabilizer[list_List,GenSet[],stableGS_]:=stableGS;

The function is not protected because it will be overloaded below with definitions on SGSs.

Using the group S3 we see the difference between pointwise stabilization and setwise stabilization:

GenSet@@Dimino[GenSet[Perm[{2,3,1}],Perm[{2,1,3}]]]

Stabilizer[{1,2},%]

SetStabilizer[{1,2},%%]

4. Strong generating sets

(*************** 4. Strong generating sets ***************)

{TimeUsed[],SessionTime[],TimeUsed[]/SessionTime[]}

4.0. Comments

The fast manipulation of permutation groups is based on the combination of two "complementary" ideas: orbits (the description of how points move under the group) and stabilization (the description of how points don't move under the group). A strong generating sets (SGS) is a clever description of a group which allows efficient encoding of both simultaneously.

From now on most algorithms are encoded using an imperative paradigm (like C), to follow Butler as closely as possible, though this is rather inefficient in Mathematica. It really doesn't matter because in the end we will recode the whole thing in C, much faster, where the imperative paradigm is natural. Therefore this Mathematica code can be considered as a testbed for new ideas and algorithms. Once they work here then they are transferred to the C code.

If[$ReadingVerbose,Print["Reading section 4: Strong generating sets"],Null,Null]

4.1. Orbits and Schreier vectors

4.1.0. Comments

An orbit is an irreducible set of points closed under OnPoints. Given a group of permutations and a point of the orbit this is a well defined concept, independent of the notation used for the permutations. However, when we want to find all orbits we need to know in advance the length of the set of points, because we want to partition it in orbits. In the same way, in order to give a Schreier vector we need to know the length of the set to construct vectors of that length. That is why Orbit does not have a length argument, but Orbits, SchreierOrbit and SchreierOrbits do have it.

4.1.1. Orbit

Orbit and Orbits give the orbits without their Schreier vectors.When no generators are given it means that we are calculating the orbit of a point under the group {ID}. Points in the orbits are not sorted and the initial point is used as first point.

This is algorithm 7.1 of Butler:

?Orbit

Orbit is simply a driver for the Mathematica or the C codes:

(* Driver *)
Options[Orbit]={MathLink:>$xpermQ};
Orbit[point_Integer?Positive,StrongGenSet[base_,GS_],options___]:=Orbit[point,GS,options];
Orbit[point_Integer?Positive,GenSet[],options___]:={point};
Orbit[point_Integer?Positive,GS_GenSet,options:OptionsPattern[]]:=If[OptionValue[MathLink],MathLinkOrbit,MathOrbit][point,GS];
Orbit[point_Integer?Positive,GS_]:=Throw@Message[Orbit::unknown,"generating set",GS];
Orbit[point_,GS_]:=Throw@Message[Orbit::unknown,"point",point];
SetNumberOfArguments[Orbit,{2,Infinity}];
Protect[Orbit];
(* Mathematica code *)
MathOrbit[point_,GS_]:=Module[{orbit={point},np=1,gamma,newgamma},
While[np≤Length[orbit],
gamma=orbit[[np]];
Map[If[FreeQ[orbit,newgamma=OnPoints[gamma,#]],AppendTo[orbit,newgamma]]&,GS];
++np];
orbit
];

Examples:

$xpermQ=False

test[Orbit,
Orbit[1,GenSet[Perm[{2,3,4,1}]]],
{1,4,3,2}
]

test[Orbit,
Orbit[5,GenSet[Rules[1->2,2->1],Rules[4->5,5->4]]],
{5,4}
]

4.1.2. Orbits

In order to give all orbits we need to give the length of the set of points to partition:

This is algorithm 7.2 of Butler:

?Orbits

(* Compute PermLength if not given *)
Orbits[gs_,initpoints_List:{}]:=Orbits[gs,PermLength[gs],initpoints];
(* Transform SGS to GS *)
Orbits[StrongGenSet[base_,GS_],len_Integer,initpoints_List:{}]:=Orbits[GS,len,initpoints];
(* Compute orbits using Orbit *)
Orbits[GS_GenSet,len_Integer,initpoints_List:{}]:=Fold[If[MemberQ[Flatten[#1],#2],#1,Append[#1,Orbit[#2,GS]]]&,{},initpoints~Join~Complement[Range[len],initpoints]];
(* Special cases *)
Orbits[(Symmetric|Antisymmetric)[list_List],len_Integer]:=Join[{list},List/@Complement[Range[len],list]];
Orbits[(Symmetric|Antisymmetric)[list_List],len_Integer,initpoints_List]:=Join[{SortB[list,initpoints]},List/@Complement[Range[len],list]];
SetNumberOfArguments[Orbits,{1,3}];
Protect[Orbits];

Recall that points in the orbits are not reordered after computation (we shall do that, however, when computing SchreierOrbit).

Examples:

test[Orbits,
Orbits[GenSet[Rules[1->2,2->1],Rules[4->5,5->4]],6],
{{1,2},{3},{4,5},{6}}
]

test[Orbits,
Orbits[GenSet[Rules[1->2,2->1],Rules[4->5,5->4]]],
{{1,2},{3},{4,5}}
]

test[Orbits,
Orbits[GenSet[],5],
{{1},{2},{3},{4},{5}}
]

test[Orbits,
Orbits[GenSet[Rules[1->2,2->1],Rules[4->5,5->4]],6,{2,5}],
{{2,1},{5,4},{3},{6}}
]

test[Orbits,
Orbits[Symmetric[{1,3,5}],6],
{{1,3,5},{2},{4},{6}}
]

test[Orbits,
Orbits[Antisymmetric[{1,3,5}],6,{5,3}],
{{5,3,1},{2},{4},{6}}
]

If the given length is too low then we can miss some of the orbits (we are asking to partition {1,2,3} only):

test[Orbits,
Orbits[GenSet[Rules[1->2,2->1],Rules[4->5,5->4]],3,{2}],
{{2,1},{3}}
]

4.1.3. SchreierOrbit

We need to give the lenght of nu, w. When we supply a SGS instead of just the GS, we follow Canon code in this convention: use as first point of the orbit the mininum of the orbit with respect to the base of the SGS. We haven´t done that in Orbit and Orbits.

The output is a Schreier structure, having three elements: the orbit, the nu vector of permutations and the w vector of backward pointers.

?Schreier

SetNumberOfArguments[Schreier,{2,Infinity}];
Protect[Schreier];

Schreier orbit:

?SchreierOrbit

Unfortunately, this line introduces a redundant computation of the orbit, but I don't see how to avoid it (MinB is, however, very efficient):

SchreierOrbit[point_Integer,StrongGenSet[B_List,GS_GenSet],deg_Integer]:=SchreierOrbit[MinB[Orbit[point,GS],B],GS,deg];

For special cases perhaps there is a faster solution than this:

SchreierOrbit[point_Integer,sym:(_Symmetric|_Antisymmetric),deg_Integer]:=SchreierOrbit[point,SchreierSims[{},sym],deg];

This is algorithm 7.4 of Butler:

SchreierOrbit::infty="Found Infinity as a point, with generating set `1`.";

SchreierOrbit[point_Integer?Positive,gs_]:=SchreierOrbit[point,gs,PermLength[gs]];
SchreierOrbit[point_Integer?Positive,GenSet[],len_Integer]:=Schreier[{point},Table[0,{len}],Table[0,{len}]];
SchreierOrbit[point_Integer?Positive,GS:GenSet[__],len_Integer]:=Module[{orbit={point},nu,w,np=1,gamma,newgamma},
nu=w=Table[0,{len}];
While[np≤Length[orbit],
gamma=orbit[[np]];
Map[If[FreeQ[orbit,newgamma=OnPoints[gamma,#]],AppendTo[orbit,newgamma];nu[[newgamma]]=#;w[[newgamma]]=gamma]&,GS];
++np];
Schreier[orbit,nu,w]
];
SchreierOrbit[Infinity,GS_,len_]:=Throw@Message[SchreierOrbit::infty,GS];
SetNumberOfArguments[SchreierOrbit,{2,3}];
Protect[SchreierOrbit];

Examples:

test[SchreierOrbit,
SchreierOrbit[1,GenSet[],4],
Schreier[{1},{0,0,0,0},{0,0,0,0}]
]

test[SchreierOrbit,
SchreierOrbit[1,GenSet[Perm[{2,3,4,1}]]],
Schreier[{1,4,3,2},{0,Perm[{2,3,4,1}],Perm[{2,3,4,1}],Perm[{2,3,4,1}]},{0,3,4,1}]
]

test[SchreierOrbit,
SchreierOrbit[5,GenSet[Rules[1->2,2->1],Rules[4->5,5->4]],6],
Schreier[{5,4},{0,0,0,Rules[4->5,5->4],0,0},{0,0,0,5,0,0}]
]

test[SchreierOrbit,
SchreierOrbit[3,StrongGenSet[{3,1},GenSet[-Cycles[{1,2}],-Cycles[{3,4}],Cycles[{1,3},{2,4}]]],4],
Schreier[{3,4,1,2},{Cycles[{1,3},{2,4}],Cycles[{1,3},{2,4}],0,-Cycles[{3,4}]},{3,4,0,3}]
]

The key point is the fact that g=OnPoints[ w[[g]], nu[[g]] ] for any g. That is, we know that we can go from w[[g]] to g using nu[[g]].

test[SchreierOrbit,
SchreierOrbit[2,StrongGenSet[{1,2,3},GenSet[Cycles[{1,2}],Cycles[{2,3}]]],3],
Schreier[{1,2,3},{0,Cycles[{1,2}],Cycles[{2,3}]},{0,1,2}]
]

4.1.4. SchreierOrbits

We can have all orbits together in a single vector:

?SchreierOrbits

(* Compute len from gs using PermLength *)
SchreierOrbits[gs_,initpoints_List:{}]:=SchreierOrbits[gs,PermLength[gs],initpoints];
(* GS: Compute all Schreier orbits from SchreierOrbit and Orbits *)
SchreierOrbits[GS:GenSet[__],len_Integer,initpoints_List:{}]:=Schreier[Sequence@@#[[1]],Plus@@#[[2]],Plus@@#[[3]]]&@Transpose[SchreierOrbit[#,GS,len]&/@First/@Orbits[GS,len,initpoints]/.Schreier->List];
(* SGS: Use base as initpoints if initpoints not specified *)
SchreierOrbits[StrongGenSet[B_List,GS:GenSet[__]],len_Integer]:=SchreierOrbits[GS,len,B];
SchreierOrbits[StrongGenSet[B_List,GS:GenSet[__]],len_Integer,initpoints_List]:=SchreierOrbits[GS,len,initpoints];
(* Special cases *)
SchreierOrbits[sym:(_Symmetric|_Antisymmetric),len_Integer,initpoints_:{}]:=SchreierOrbits[GenSet[sym],len,initpoints];
SchreierOrbits[GenSet[],len_Integer,initpoints_:{}]:=Schreier[Sequence@@Orbits[GenSet[],len],Table[0,{len}],Table[0,{len}]];
SchreierOrbits[StrongGenSet[_,GenSet[]],len_Integer]:=Schreier[Sequence@@Orbits[GenSet[],len],Table[0,{len}],Table[0,{len}]];
SetNumberOfArguments[SchreierOrbits,{1,3}];
Protect[SchreierOrbits];

Examples. Note that the number of arguments of Schreier is finally #orbits + 2.

TracePrint[SchreierOrbits[GenSet[]],SchreierOrbits[___]]

test[SchreierOrbits,
SchreierOrbits[GenSet[],5],
Schreier[{1},{2},{3},{4},{5},{0,0,0,0,0},{0,0,0,0,0}]
]

test[SchreierOrbits,
SchreierOrbits[GenSet[Rules[1->2,2->1],Rules[4->5,5->4]]],
Schreier[{1,2},{3},{4,5},{0,Rules[1->2,2->1],0,0,Rules[4->5,5->4]},{0,1,0,0,4}]
]

test[SchreierOrbits,
SchreierOrbits[GenSet[Rules[1->2,2->1],Rules[4->5,5->4]],5,{2,5}],
Schreier[{2,1},{5,4},{3},{Rules[1->2,2->1],0,0,Rules[4->5,5->4],0},{2,0,0,5,0}]
]

test[SchreierOrbits,
SchreierOrbits[StrongGenSet[{1,4},GenSet[Rules[1->2,2->1],Rules[4->5,5->4]]],5],
Schreier[{1,2},{4,5},{3},{0,Rules[1->2,2->1],0,0,Rules[4->5,5->4]},{0,1,0,0,4}]
]

test[SchreierOrbits,
SchreierOrbits[StrongGenSet[{2,5},GenSet[Rules[1->2,2->1],Rules[4->5,5->4]]],5],
Schreier[{2,1},{5,4},{3},{Rules[1->2,2->1],0,0,Rules[4->5,5->4],0},{2,0,0,5,0}]
]

Example from page 59 of Butler's book:

a=Cycles[{1,4,5,11,6,10,3,2},{7,8}];
b=Cycles[{1,4,5,11,6,10,3,2},{8,9}];
c=Cycles[{1,4,2,3},{5,11,10,6}];

test[Orbit,
Orbit[7,GenSet[a,b,c]],
{7,8,9}
]

test[SchreierOrbit,
SchreierOrbit[7,GenSet[a,b,c],11]/.{a->"a",b->"b",c->"c"},
Schreier[{7,8,9},{0,0,0,0,0,0,0,"a","b",0,0},{0,0,0,0,0,0,0,7,8,0,0}]
]

test[Orbit,
Orbit[1,GenSet[a,b,c]],
{1,4,5,2,11,3,6,10}
]

test[SchreierOrbit,
SchreierOrbit[1,GenSet[a,b,c]]/.{a->"a",b->"b",c->"c"},
Schreier[{1,4,5,2,11,3,6,10},{0,"c","c","a","a","a",0,0,0,"c","a"},{0,4,2,1,4,11,0,0,0,11,5}]
]

test[Orbits,
Orbits[GenSet[a,b,c]],
{{1,4,5,2,11,3,6,10},{7,8,9}}
]

test[SchreierOrbits,
SchreierOrbits[GenSet[a,b,c],11]/.{a->"a",b->"b",c->"c"},
Schreier[{1,4,5,2,11,3,6,10},{7,8,9},{0,"c","c","a","a","a",0,"a","b","c","a"},{0,4,2,1,4,11,0,7,8,11,5}]
]

test[SchreierOrbits,
SchreierOrbits[GenSet[a,b,c],11,{2,8}]/.{a->"a",b->"b",c->"c"},
Schreier[{2,1,3,4,5,11,6,10},{8,7,9},{"a",0,"c","a","a","a","a",0,"b","c","a"},{2,0,2,1,4,11,8,0,8,11,5}]
]

If no permutations are given, we just return the trivial sets. Note the relevance of the len argument:

test[Orbit,
Orbit[1,GenSet[]],
{1}
]

test[SchreierOrbit,
SchreierOrbit[1,GenSet[],3],
Schreier[{1},{0,0,0},{0,0,0}]
]

test[Orbits,
Orbits[GenSet[],3],
{{1},{2},{3}}
]

test[SchreierOrbits,
SchreierOrbits[GenSet[],3],
Schreier[{1},{2},{3},{0,0,0},{0,0,0}]
]

4.1.5. TraceSchreier

Finally, we can trace the Schreier vectors, which is the operation giving the permutation moving the first-point of an orbit to a given point:

?TraceSchreier

This is algorithm 7.3 of Butler:

TraceSchreier::trace="Cannot trace `1` with orbit `2`";
TraceSchreier[p_Integer,oS:Schreier[orbits___List,nu_,w_]]:=Which[
MemberQ[First/@{orbits},p],ID,
p===0,Throw[Message[TraceSchreier::trace,0,oS]],
Length[nu]<p,ID,
True,PermProduct[TraceSchreier[w[[p]],oS],nu[[p]]]
];
TraceSchreier[p_Integer,os_Schreier]:=Throw@Message[TraceSchreier::trace,p,os];
SetNumberOfArguments[TraceSchreier,2];
Protect[TraceSchreier];

Examples:

test[SchreierOrbit,
SchreierOrbit[1,GenSet [a,b,c],11]/.{a->"a",b->"b",c->"c"},
Schreier[{1,4,5,2,11,3,6,10},{0,"c","c","a","a","a",0,0,0,"c","a"},{0,4,2,1,4,11,0,0,0,11,5}]
]

test[TraceSchreier,
TraceSchreier[11,SchreierOrbit[1,GenSet [a,b,c],11]],
Cycles[{1,11,3,4,6,2,5,10},{7,8}]
]

test[PermPower,
PermPower[a,3],
Cycles[{1,11,3,4,6,2,5,10},{7,8}]
]

test[TraceSchreier,
TraceSchreier[11,SchreierOrbits[GenSet[a,b,c],11]],
Cycles[{1,11,3,4,6,2,5,10},{7,8}]
]

test[OnPoints,
OnPoints[1,Cycles[{1,11,3,4,6,2,5,10},{7,8}]],
11
]

test[TraceSchreier,
TraceSchreier[1,SchreierOrbit[1,GenSet[a,b,c],11]],
ID
]

Catch@TraceSchreier[8,SchreierOrbit[1,GenSet[a,b,c],11]]

Note the following particular cases:

test[TraceSchreier,
TraceSchreier[1,Schreier[{1},{0,0,0},{0,0,0}]],
ID
]

Catch@TraceSchreier[2,Schreier[{1},{0,0,0},{0,0,0}]]

4.2. StabilizerChain

4.2.0. Comments

In the MathLink connection section we shall give a definition for StrongGenSet with three arguments.

SetNumberOfArguments[StrongGenSet,{2,3}];
Protect[StrongGenSet];

Now we deal with the stabilizing part of the SGSs. The key idea is that a SGS allows the construction of a hierarchy of subgroups: the stabilizers of consecutive points in the base.

4.2.1. StabilizerChain

This function gives the chain of stabilizers associated to a SGS. There is no check of validity of the initial pair (base, GS):

?StabilizerChain

StabilizerChain[StrongGenSet[base_List,GS_GenSet]]:=StrongGenSet[DeleteCases[base,Alternatives@@#],Stabilizer[#,GS]]&/@FoldList[Append[#1,#2]&,{},base];
StabilizerChain[sym:(_Symmetric|_Antisymmetric)]:=StabilizerChain[SchreierSims[{},sym]];
SetNumberOfArguments[StabilizerChain,1];
Protect[StabilizerChain];

If there are no generators, the orbits have trivial structure:

StabilizerChain[StrongGenSet[{1,2},GenSet[]]]//ColumnForm

StabilizerChain[StrongGenSet[{},GenSet[]]]//ColumnForm

Note that the first SGS is the original SGS and the last one is the empty SGS (generating Group[ID]):

StabilizerChain[StrongGenSet[{1,3},GenSet[Perm[{2,1,3,4,6,5}],Perm[{1,2,4,3,6,5}],Perm[{3,4,1,2,5,6}]]]]//ColumnForm

StabilizerChain[StrongGenSet[{1,2,3,4},GenSet[Perm[{2,1,3,4,6,5}],Perm[{1,2,4,3,6,5}],Perm[{3,4,1,2,5,6}]]]]//ColumnForm

4.3. Inductive foundation

4.3.0. Comments

Here we construct a number of commands which perform simple tasks using strong generating sets.

This is essentially chapter 10 of Butler's book. There are lots of examples there, which have been repeated (and some of them corrected) in the companion notebook ButlerExamples.nb.

Recall that given the SGS of a group G with base {b1, b2, ..., bk} we can construct a hierarchy of stabilizers
			{id} = Gk+1 < Gk < ... < G2 < G1 = G
where Gi is the stabilizer of the points {b1, ..., bi-1}. The notation is horrible because it should have been G0 = G and Gk = {id}, etc.

The group G1=G is divided in the subroup G2 (the stabilizer of b1) and its cosets. There is a one-to-one correspondence between those cosets and the points of the orbit b1^G. We select the smallest permutation u of those cosets, that will represent the coset in all of the following computations. Given any permutation g in G we compute b1^g; if we get b1 then that means that g is in G2 (the stabilizer of b1), but if we get another point in the orbit of b1 then we know that the permutation g belongs to the corresponding right coset G2 x g. If u is the representant of that coset, then h = g x u^-1 belongs to G2. Note that b1^g = b1^(h x u) = (b1^h)^u = b1^u. Now we can repeat the same process starting with the permutation h and the group G2, for which we have a SGS with base {b2, b3, ..., bk}. The length of the orbit b1^G is the number of cosets G:G2. Then the length of the orbit b2^G2 is the number of cosets G2:G3, etc. Finally the length of the orbit bk^Gk is the number of cosets Gk:{id}, which is the order of Gk. The orbit Delta^(i) = bi^Gi is called the basic orbit i and its length is called the basic index i.

4.3.1. OrderOfGroup

This recursive function just multiplies the lengths of the basic orbits of the SGSs in the chain of stabilizers:

?OrderOfGroup

Note that there are no options in OrderOfGroup, but we can pass the MathLink option to Orbit:

OrderOfGroup::nobase="Found strong generating set with empty base but generators `1`. Assigning ficticious order 0.";
OrderOfGroup[StrongGenSet[{___},GenSet[]],options___]:=1;
OrderOfGroup[StrongGenSet[{},GS_GenSet]]:=(Message[OrderOfGroup::nobase,GS];0);
OrderOfGroup[StrongGenSet[base_List,GS_GenSet],options___]:=Length[Orbit[First[base],GS,options]]OrderOfGroup[StrongGenSet[Rest[base],Stabilizer[{First[base]},GS]],options];
OrderOfGroup[(Symmetric|Antisymmetric)[list_List]]:=Length[list]!;
SetNumberOfArguments[OrderOfGroup,{1,Infinity}];
Protect[OrderOfGroup];

Examples:

test[OrderOfGroup,
OrderOfGroup[StrongGenSet[{1,3},GenSet[Perm[{2,1,3,4,6,5}],Perm[{1,2,4,3,6,5}],Perm[{3,4,1,2,5,6}]]]],
8
]

If there are no generators the group is just {ID}:

test[OrderOfGroup,
OrderOfGroup[StrongGenSet[{},GenSet[]]],
1
]

Note that redundant points in the base do not change the result:

test[OrderOfGroup,
OrderOfGroup[StrongGenSet[{1,2,3,4},GenSet[Perm[{2,1,3,4,6,5}],Perm[{1,2,4,3,6,5}],Perm[{3,4,1,2,5,6}]]]],
8
]

test[OrderOfGroup,
OrderOfGroup[StrongGenSet[{1,2},GenSet[]]],
1
]

These are the SGSs involved in the canonicalization of products of n (anti)symmetric 2-tensors:

sgs[n_]:=StrongGenSet[Range[1,2n-1,2],GenSet@@Join[Table[Cycles[{i,i+1}],{i,1,2n-1,2}],Table[Cycles[{i,i+2},{i+1,i+3}],{i,1,2n-3,2}]]]

sgs[3]

OrderOfGroup[sgs[#]]&/@Range[25]

%//N

4.3.2. Membership. PermMemberQ, PermWord

Recursive function that checks membership of a permutation in a group generated by a SGS:

This is algorithm 10.2 of Butler. Algorithms 10.1 and 10.3 (which produce the same output) have not been implemented in xPerm.

?PermMemberQ

This function is iterative, with output {residue, word}, and needs to store the previous stages of the computation, in the third argument.

?PermWord

PermWord[perm_,StrongGenSet[base_List,GS_GenSet],word_List:{}]:=If[Length[base]===0,
Prepend[word,perm],
Module[{
Sorbit=SchreierOrbit[First[base],GS,Max[PermDeg[perm],PermDeg[GS]]],
point=OnPoints[First[base],perm],
u},
If[MemberQ[First[Sorbit],point],
u=TraceSchreier[point,Sorbit];
PermWord[PermProduct[perm,InversePerm[u]],StrongGenSet[Rest[base],Stabilizer[{First[base]},GS]],Prepend[word,u]],
Prepend[word,perm]]
]
];
PermMemberQ[perm_,sym:(_Symmetric|_Antisymmetric)]:=PermMemberQ[perm,SchreierSims[{},sym]];
PermMemberQ[perm_,SGS_]:=With[{id=ID[perm],first=First@PermWord[perm,SGS]},Or[PermEqual[first,id],PermEqual[first===-id]&&MemberQ[SGS[[2]],first]]];
SetNumberOfArguments[PermWord,{2,3}];
SetNumberOfArguments[PermMemberQ,2];
Protect[PermWord,PermMemberQ];

Special cases:

test[PermMemberQ,
PermMemberQ[ID,StrongGenSet[{1},GenSet[]]],
True
]

test[PermMemberQ,
PermMemberQ[ID,StrongGenSet[{},GenSet[]]],
True
]

test[PermMemberQ,
PermMemberQ[Perm[{2,1}],StrongGenSet[{1},GenSet[]]],
False
]

test[PermMemberQ,
PermMemberQ[Perm[{2,1}],StrongGenSet[{},GenSet[]]],
False
]

test[PermMemberQ,
PermMemberQ[Rules[],StrongGenSet[{},GenSet[]]],
True
]

test[PermWord,
PermWord[ID,StrongGenSet[{1},GenSet[]]],
{ID,ID}
]

test[PermWord,
PermWord[ID,StrongGenSet[{},GenSet[]]],
{ID}
]

test[PermWord,
PermWord[Perm[{2,1}],StrongGenSet[{1},GenSet[]]],
{Perm[{2,1}]}
]

test[PermWord,
PermWord[Perm[{2,1}],StrongGenSet[{},GenSet[]]],
{Perm[{2,1}]}
]

test[PermWord,
PermWord[Rules[],StrongGenSet[{},GenSet[]]],
{Rules[]}
]

4.3.3. FromBaseImage

We know that a permutation is uniquely defined by the images of a base under that permutation. This function finds the permutation corresponding to a given list of images and a SGS.

This is algorithm 10.5 of Butler.

?FromBaseImage

FromBaseImage::noimgs="Invalid list of images `1`.";
FromBaseImage[images_List,sym:(_Symmetric|_Antisymmetric),len___]:=FromBaseImage[images,SchreierSims[{},sym],len];
FromBaseImage[images_List,SGS_StrongGenSet]:=FromBaseImage[images,SGS,PermLength[SGS]];
FromBaseImage[{},StrongGenSet[{},GenSet[]],0]:=ID;
FromBaseImage[images_List,SGS:StrongGenSet[base_List,GS_GenSet],len_Integer]:=Module[{i,u,g=ID,schvecs=SchreierOrbits[#,len]&/@StabilizerChain[SGS],imgs=images},
If[Length[images]=!=Length[base],Throw[Message[FromBaseImage::noimgs,images]]];
For[i=1,i≤Length[imgs],i++,
(* If[MemberQ[schvecs[[i,1]],imgs[[i]]],*)
u=TraceSchreier[imgs[[i]],schvecs[[i]]];
g=PermProduct[u,g];
imgs=OnPoints[imgs,InversePerm[u]]
(* ,
Throw[Message[FromBaseImage::noimgs,images]];
]*)
];
If[imgs==base,g,Throw[Message[FromBaseImage::noimgs,images]]]];
SetNumberOfArguments[FromBaseImage,{2,3}];
Protect[FromBaseImage];

Example:

GS=GenSet[-Cycles[{1,2}],-Cycles[{3,4}],-Cycles[{5,6}],-Cycles[{7,8}],-Cycles[{9,10}],Cycles[{1,3},{2,4}],Cycles[{3,5},{4,6}],Cycles[{5,7},{6,8}],Cycles[{7,9},{8,10}]];
base={1,3,5,7,9};
SGS=StrongGenSet[base,GS];

perm=PermProduct@@GS

OnPoints[base,perm]

test[FromBaseImage,
FromBaseImage[%,SGS],
-Cycles[{1,10,7,6,3,2,9,8,5,4}]
]

word=PermWord[-Cycles[{1,10,7,6,3,2,9,8,5,4}]
,SGS]

PermProduct@@word

%===perm

Let us have a closer look at that example. The SGS gives the following chain of stabilizers:

(chain=StabilizerChain[SGS])//ColumnForm

OrderOfGroup/@chain

The last permutation in the word is a permutation of the group, but not the first stabilizer. Previous perms in the word belong to more stabilizers. The last line shows simply the fact that the first element of the word (ID) belongs to the last stabilizer (the identity group).

Outer[PermMemberQ,Reverse[word],chain]//MatrixForm

If the list of images does not correspond to an element of the group, then we get an error message:

Catch@FromBaseImage[{5,2,3,9,4},StrongGenSet[base,GS],10]

We can get a residue:

perm=Cycles[{1,2,3,4,5,6,7}];

PermWord[perm,SGS]

PermProduct@@%

such that next permutation does belong to the group:

PermProduct@@Rest[%%]

PermMemberQ[%,SGS]

perm=.
GS=.
base=.
SGS=.
word=.

4.3.4. Enumeration of all base images

This is algorithm 10.8 of Butler.

?AllBaseImages

This function lists all elements of a group. Note that we use PermDeg and not PermLength because it produces shorter or equal Schreier vectors.

AllBaseImages[sym:(_Symmetric|_Antisymmetric)]:=AllBaseImages[SchreierSims[{},sym]];
AllBaseImages[SGS:StrongGenSet[base_List,GS_GenSet]]:=AllBaseImages[SGS,PermDeg[SGS]];
AllBaseImages[SGS:StrongGenSet[base_List,GS_GenSet],len_Integer]:=
BaseImage[1,Length[base],{},{ID},SchreierOrbits[#,len]&/@StabilizerChain[SGS]];
SetNumberOfArguments[AllBaseImages,{1,2}];
Protect[AllBaseImages];

Internal recursive part:

BaseImage[i_,k_,points_List,word_List,schvecs_]:=Module[{g=PermProduct@@word,ig,j,u,Deltab},
If[i==k+1,{Rule[points,g]},
ig=InversePerm[g];
Deltab=OnPoints[schvecs[[i,1]],g];
Flatten[Table[u=TraceSchreier[OnPoints[Deltab[[j]],ig],schvecs[[i]]];
BaseImage[i+1,k,Append[points,Deltab[[j]]],{u,g},schvecs],
{j,1,Length[Deltab]}]]
]
];

Examples. Note that not all combinations of images are allowed. For instance images {1,2} do not identify a valid permutation:

AllBaseImages[StrongGenSet[{1,3},GenSet[-Cycles[{1,2}],-Cycles[{3,4}],Cycles[{1,3},{2,4}]]]]//ColumnForm

AllBaseImages[Antisymmetric[{1,2,3}]]//Column

Catch@FromBaseImage[{1,2},StrongGenSet[{1,3},GenSet[-Cycles[{1,2}],-Cycles[{3,4}],Cycles[{1,3},{2,4}]]]]

Note the order of the permutations. Ordering the bases with normal ordering, as we have just done produces also normal ordering in the group:

group=Dimino[GenSet[-Cycles[{1,2}],-Cycles[{3,4}],Cycles[{1,3},{2,4}]]]

PermSort[group]

However, sorting with the base produces a different output:

PermSort[group,{1,3}]

v9 note : I have the impression that the ordering is not correct now...

Special cases:

AllBaseImages[StrongGenSet[{},GenSet[]]]//ColumnForm

AllBaseImages[StrongGenSet[{1},GenSet[]]]

The symmetric group. Here all base images are allowed:

SGS=StrongGenSet[{1,2,3},GenSet[Perm[{2,1,3,4}],Perm[{1,3,2,4}],Perm[{1,2,4,3}]]];

AllBaseImages[SGS]//ColumnForm

FromBaseImage[{4,1,3},SGS]

SGS=.

Let us now compare timings of Dimino and AllBaseImages:

SGS=StrongGenSet[{1,2,3,4,5,6},GenSet[Perm[{2,1,3,4,5,6,7}],Perm[{1,3,2,4,5,6,7}],Perm[{1,2,4,3,5,6,7}],Perm[{1,2,3,5,4,6,7}],Perm[{1,2,3,4,6,5,7}],Perm[{1,2,3,4,5,7,6}],Perm[{7,2,3,4,5,6,1}]]];

AbsoluteTiming[permsMathematica=Permutations[{1,2,3,4,5,6,7}];]

AbsoluteTiming[permsDimino=Dimino[SGS[[2]]];]

AbsoluteTiming[permsAllBases=AllBaseImages[SGS];]

{Length[permsMathematica],Length[permsDimino],Length[permsAllBases]}

4.4. Reordering points. Base-based ordering

Example:

GS=GenSet[Cycles[{1,2,3},{4,5,6},{7,8,9}],Cycles[{2,4,3,7},{5,6,9,8}],Cycles[{2,5,3,9},{4,8,7,6}],Cycles[{1,10},{4,5},{6,8},{7,9}],Cycles[{1,11},{4,6},{5,9},{7,8}]];
base={11,10,1,2};
SGS=StrongGenSet[base,GS];

OrderOfGroup[SGS]

Take two permutations and their base images:

perm1=PermProduct@@GS

baseimages1=OnPoints[base,perm1]

perm2=PermProduct@@GS[[First@RandomPerm[5]]]

baseimages2=OnPoints[base,perm2]

Now we can decide whether the permutations are sorted or not:

comp[i_,i_]:=Null;
comp[i_,j_]:=Throw[i<j];
Catch@Inner[comp,baseimages1,baseimages2,List]

This is equivalent to lexicographic order of the permutations if we reorder the points so that the beginning of the set is given by the base:

points=Join[base,Complement[Range[11],base][[First@RandomPerm[7]]]]

images1=OnPoints[points,perm1]

images2=OnPoints[points,perm2]

Catch@Inner[comp,images1,images2,List]

Note that the ordering of the points before mapping is independent of the order of points after mapping. In this example, for instance, we have used the ordering given by the base before mapping, and then the usual numeric order for integers.

However, the question is whether any ordering of the points after mapping has the identity as first element of a set of sorted permutations. We see that in our example that is false:

perm0=ID;

baseimages0=OnPoints[base,perm0]

Now we can decide whether the permutations are sorted or not:

Catch@Inner[comp,baseimages0,baseimages1,List]

This is equivalent to lexicographic order of the permutations if we reorder the points so that the beginning of the set is given by the base:

images0=OnPoints[points,perm0]

Catch@Inner[comp,images0,images1,List]

Now we change to a new ordering based on the base, and now the identity is first:

compB[i_,i_]:=Null;
compB[i_,j_]:=Throw[LessB[{i,j},base]];
Catch@Inner[compB,baseimages0,baseimages1,List]

Catch@Inner[compB,images0,images1,List]

4.5. Backtrack search

4.5.0. Comments

Given a group G with a SGS we want to construct a subgroup H of G identifying another sgs for H with the same base. The subgroup H is identified by the decidable property P. The algorithm proceeds constructing subgroups Hi of the stabilizers Gi from the smallest stabilizers Gk, Gk-1, etc towards the largest stabilizers G2, G1. We first take the whole Gk-1 and construct Hk-1, checking P on all permutations of Gk-1. Now we go to Gk-2. The key point is that we do not have to check all permutations in Gk-2 \ Gk-1 because we can use the cosets formed from the subgroup Hk-1. For example, if g is a permutation of Gk-2 obeying (or not obeying) the property P then it is clear that all permutations in Hk-1 x g will obey it too (or will not obey it either). This can be easily implemented through the use of base images.

4.5.1. Search and Generate

The group G is described by its stabilizer-chain. The procedure is divided in two parts: there is a first routine constructing each stabilizer of the group K (function Search) and a recursive routine finding new permutations of their generating sets (funtion Generate). The "property" argument is a Boolean function acting on a permutation of G.

?Search

Options[Search]={xPermVerbose->False};
Search[SGS_StrongGenSet,property_,s_Integer,SGSK_,options___]:=Search[StabilizerChain[SGS],property,s,SGSK,options];
Search[chain_List,property_,s_Integer,SGSK_,options:OptionsPattern[]]:=Module[{base=chain[[1,1]],newSGSK=SGSK,k=Length[chain]-1,Deltas,gammas,Korbit,verb},
verb=OptionValue[xPermVerbose];
If[s==k+1,
newSGSK=StrongGenSet[base,GenSet[]],
newSGSK=Search[chain,property,s+1,SGSK,options];
Deltas=Orbit[base[[s]],chain[[s,2]]];
If[verb,Print["Branching over points ",Deltas]];
(* Avoid rechecking permutations. Not in Butler's algorithm *)
If[s=!=k,Deltas=Drop[Deltas,1]];
Do[
gammas=Deltas[[j]];
Korbit=Orbit[gammas,newSGSK[[2]]];
If[gammas===MinB[Korbit,base],
newSGSK=Generate[chain,property,s,s+1,Append[Take[base,s-1],gammas],newSGSK,verb];
If[Head[newSGSK]===Times,newSGSK=-newSGSK]]
,{j,1,Length[Deltas]}];
];
newSGSK
];
SetNumberOfArguments[Search,{4,Infinity}];
Protect[Search];

When a new permutation is added to the SGS of K then we must return to Search, independently from the recursion level of Generate. To do that we use a dirty trick: marking the SGS of K with a minus sign in front.

Generate[chain_List,property_,s_,i_,list_,SGSK_,verb_]:=Module[{base=chain[[1,1]],k=Length[chain]-1,g,newSGSK=SGSK,Deltag,gammai,otherSGSK},
g=FromBaseImage[list,ReplacePart[chain[[1]],Take[base,Length[list]],1]];
If[i==k+1,
If[verb,Print["Generate at level ",s," with i=",i,". We have list ",list," and permutation ",g]];
If[Not@PermEqual[g,ID]&&property[g],
newSGSK=-StrongGenSet[First@SGSK,Append[Last[SGSK],g]];
If[verb,Print["  Added permutation ",g]]],
Deltag=OnPoints[Orbit[base[[i]],chain[[i,2]]],g];
If[verb,Print["Generating over points ",Deltag]];
Do[
gammai=Deltag[[j]];
newSGSK=Generate[chain,property,s,i+1,Append[list,gammai],newSGSK,verb];
If[Head[newSGSK]===Times,Break[]],
{j,1,Length[Deltag]}];
];
newSGSK
];

Butler has additional improvements based on restrictions on image points, but we do not implement that here.

Example: All permutations:

Search[StrongGenSet[{1,2,3,4},GenSet[Perm[{2,1,3,4,5}],Perm[{1,3,2,4,5}],Perm[{1,2,4,3,5}],Perm[{1,2,3,5,4}],Perm[{5,2,3,4,1}]]],True&,1,StrongGenSet[{},GenSet[]],xPermVerbose->True]

OrderOfGroup[%]

4.5.2. Examples: Centralizer

There are two examples given by Butler in pages 100 and 101:

z=Cycles[{1,3},{2,4}];

SGS=Search[StrongGenSet[{1,2,3},GenSet[Cycles[{1,2}],Cycles[{2,3}],Cycles[{3,4}],Cycles[{1,4}]]],PermEqual[PermProduct[z,#],PermProduct[#,z]]&,1,StrongGenSet[{},GenSet[]],xPermVerbose->True]

StabilizerChain[SGS]//ColumnForm

OrderOfGroup[SGS]

Dimino[SGS[[2]]]

%//Length

z=Cycles[{1,2},{3,4}];

SGS=Search[StrongGenSet[{1,2,3},GenSet[Cycles[{1,2}],Cycles[{2,3}],Cycles[{3,4}],Cycles[{1,4}]]],PermEqual[PermProduct[z,#],PermProduct[#,z]]&,1,StrongGenSet[{},GenSet[]],xPermVerbose->True]

StabilizerChain[SGS]//ColumnForm

OrderOfGroup[SGS]

Dimino[SGS[[2]]]

%//Length

4.6. BaseChange, DeleteRedundantGenerators

We want to return a new SGS with points j and j+1 exchanged in base.

Generalized basic orbit:

BasicOrbit[j_Integer,k_Integer,StrongGenSet[base_List,GS_GenSet]]:=Orbit[base[[j]],Stabilizer[Take[base,k],GS]]/;k<j;

This is algorithm 12.3 of Butler. Algorithms 12.1 and 12.2 are not programmed in xPerm. It is private:

Interchange[SGS:StrongGenSet[base_List,GS_GenSet],j_Integer,len_Integer]:=Module[{Deltaj,Deltajp1,LDeltaBarjp1,T,Gamma,Delta,gamma,p,g1,g2},
Deltaj=BasicOrbit[j,j-1,SGS];
Deltajp1=BasicOrbit[j+1,j,SGS];
LDeltaBarjp1=Length[Deltajp1]Length[Deltaj]/Length[BasicOrbit[j+1,j-1,SGS]];
T=Stabilizer[Take[base,j+1],GS];
Gamma=Complement[Deltaj,base[[{j,j+1}]]];
Delta={base[[j]]};
While[Length[Delta]≠LDeltaBarjp1,
gamma=First[Gamma];
g1=TraceSchreier[gamma,SchreierOrbit[base[[j]],StrongGenSet[base[[Range[j,Length[base]]]],Stabilizer[base[[Range[j-1]]],GS]],len]];
p=OnPoints[base[[j+1]],InversePerm[g1]];
If[MemberQ[Deltajp1,p],
g2=TraceSchreier[p,SchreierOrbit[base[[j+1]],StrongGenSet[base[[Range[j+1,Length[base]]]],Stabilizer[base[[Range[j]]],GS]],len]];
AppendTo[T,PermProduct[g2,g1]];
Delta=Orbit[base[[j]],T];
Gamma=Complement[Gamma,Delta],
Gamma=Complement[Gamma,Orbit[gamma,T]]
]
];
StrongGenSet[base/.{base[[j]]->base[[j+1]],base[[j+1]]->base[[j]]},Union[GS,T]]
]

Example (page 121 of Butler):

base={1,9,8,10,2,12};
GS=GenSet[Cycles[{1,8,9},{2,11,15},{3,10,12},{4,14,19},{5,16,17},{6,21,20},{7,13,18}],Cycles[{9,18,20},{12,19,17}],Cycles[{10,21,11},{13,16,14}],Cycles[{8,13,21},{10,14,16}],Cycles[{2,6,3},{4,5,7}],Cycles[{12,20,15},{17,19,18}]];
SGS=StrongGenSet[base,GS];

There are six elements in the base:

BasicOrbit[#,#-1,SGS]&/@Range[6]

Length/@%

Times@@%

OrderOfGroup[SGS]

newSGS=Interchange[SGS,1,21]

Complement[Last[newSGS],GS]

newSGS=Interchange[newSGS,2,21]

OrderOfGroup[newSGS]

Complement[Last[newSGS],GS]

This is the complete algorithm. This is 12.5 of Butler. Note that the new base is not necessarily a simple reordering of the old base. It can have new points.

?BaseChange

BaseChange[SGS_,base_]:=BaseChange[SGS,base,PermLength[SGS]];
BaseChange[StrongGenSet[{},GenSet[]],base_,len_]:=StrongGenSet[base,GenSet[]];
BaseChange[SGS:StrongGenSet[{b_,___},_],{b_},len_]:=SGS;
BaseChange[SGS:StrongGenSet[base_List,GS_GenSet],newbase_List,len_]:=Module[{i=0,j,g=ID,more,B=base,gs=GS,pos,newSGS},
more=If[Length[newbase]>0,
MemberQ[BasicOrbit[1,0,SGS],newbase[[1]]],
False];
While[more,
i=i+1;
g=PermProduct[TraceSchreier[OnPoints[newbase[[i]],InversePerm[g]],SchreierOrbits[SGS,len]],g];
more=(i+1)≤Min[Length[B],Length[newbase]];
If[more,more=MemberQ[OnPoints[newbase[[i+1]],InversePerm[g]],BasicOrbit[i+1,i,SGS]]]
];
If[Not@PermEqual[g,ID],B=OnPoints[B,g];gs=Map[PermProduct[InversePerm[g],#,g]&,gs]];
For[j=i+1,j≤Length[newbase],j++,
If[MemberQ[B,newbase[[j]]],pos=Position[B,newbase[[j]]][[1,1]],AppendTo[B,newbase[[j]]];pos=Length[B]];
While[pos≠j &&pos>1,newSGS=Interchange[StrongGenSet[B,gs],pos-1,len];B=newSGS[[1]];gs=newSGS[[2]];pos=pos-1];
];
StrongGenSet[B,gs]
];
SetNumberOfArguments[BaseChange,{2,3}];
Protect[BaseChange];

BaseChange[SGS,{1},21]//Trace

BaseChange[SGS,{9,8},21]

OrderOfGroup[%]

The solution given by Butler is wrong because he gives a permutation which is not in the original group:

PermMemberQ[Cycles[{1,7,6},{4,5,11}],SGS]

This is essentially algorithm 12.4 of Butler.  It is more complicated because Butler's algorithm's efficiency depends on the order of the generators in GS. We introduce a Permutations structure to consider all possible orderings of generators and take the most convenient case.

?DeleteRedundantGenerators

DeleteRedundantGenerators[StrongGenSet[base_,GS_GenSet]]:=Module[{Si,Sip1,i,T=GenSet[],orbit,gens,pgens,check,toadd},
Sip1=GenSet[];
check[i_,gen_,t_]:=FreeQ[Orbit[base[[i]],t],OnPoints[base[[i]],gen]];
toadd[i_,gs_,T_]:=Module[{TT=T},If[check[i,#,TT],AppendTo[TT,#]]&/@gs;Complement[TT,T]];
For[i=Length[base],i≥1,i--,
Si=Stabilizer[Take[base,i-1],GS];
gens=Complement[Si,Sip1];
pgens=Permutations[gens];
T=Join[T,First@Sort[toadd[i,#,T]&/@pgens]];
Sip1=Si;
];
StrongGenSet[base,T]
];
SetNumberOfArguments[DeleteRedundantGenerators,1];
Protect[DeleteRedundantGenerators];

Examples:

Complement[Last[newSGS],Last[DeleteRedundantGenerators[newSGS]]]

base=.
GS=.
SGS=.
newSGS=.

4.7. Schreier-Sims methods 3 and 4

Notations changed from Butler's book:
	current_gens : Si
	old_gens : oldSi

Assuming that S^(i) is a GS for H^(i) and that B,S^(i+1) are base and SGS for H^(i+1), find base and SGS for H^(i). T is the set of additional generators in S^(i) since the previous call to the procedure with the present value of i. Assume that a base and SGS set of <S^(i)-T>, (the previous value of H^(i)) are included in B and S. The present value of nu^(i) must be an extension of the previous value. This function is private.

SchreierSimsStep[B_List,S_GenSet,len_Integer,i_Integer,T_GenSet,options:OptionsPattern[]]:=Module[{Si,oldSi,genset,Deltai,oldDeltai,gn,j,k,gamma,sn,s,g,gbar,newB,newS,level,n,rules,method,verb,word,junk},

{rules,method,verb}=OptionValue[SchreierSims,{options},{UseRules,Method,xPermVerbose}];

newB=B;newS=S;
If[verb,Print["Schreier-Sims-Step called with B=",B,", S=",S/.rules,", i=",i,", T=",T/.rules]];

(* Original generating sets *)
Si=Stabilizer[Take[B,i-1],S];
oldSi=Complement[Si,T];
(* Basic orbits *)
oldDeltai=SchreierOrbit[B[[i]],oldSi,len];
Deltai=SchreierOrbit[B[[i]],Si,len];
(* Check that Deltai is an extension of oldDeltai *)
For[n=1,n≤Length[oldDeltai[[2]]],++n,
If[(Deltai[[2,n]]=!=oldDeltai[[2,n]])&&(oldDeltai[[2,n]]=!=0),
If[verb,Print["Modifying orbit Delta to extend previous value."]];
Deltai[[2,n]]=oldDeltai[[2,n]];
Deltai[[3,n]]=oldDeltai[[3,n]]
]
];

(* Loop over elements gamma of basic orbit *)
For[gn=1,gn≤Length[Deltai[[1]]],++gn,gamma=Deltai[[1,gn]];

If[MemberQ[First[oldDeltai],gamma],genset=T,genset=Si];

(* Loop sn over generators s in genset *)
For[sn=1,sn≤Length[genset],++sn,s=genset[[sn]];

(* Compute Schreier generator *)
g=PermProduct[TraceSchreier[gamma,Deltai],s,InversePerm[TraceSchreier[OnPoints[gamma,s],Deltai]]];++num;

(* If g is not in subgroup H^(i+1) *)
If[Not@PermMemberQ[g,StrongGenSet[Drop[newB,i],Stabilizer[Take[newB,i],newS]]],
k=Length[newB];
Switch[method,
"Butler3",
gbar=g;
j=1;While[j≤k&&OnPoints[newB[[j]],gbar]===newB[[j]],j++],
"Butler4",
word=PermWord[g,StrongGenSet[Drop[newB,i],Stabilizer[Take[newB,i],newS]]];
gbar=First[word];
j=i+Length[word],
_,
Throw @Print["Unknown Schreier-Sims variant method: ",method]
];
(* Extend GS *)
AppendTo[newS,gbar];
If[verb,Print["Appended S[[",Length[newS],"]]=",gbar/.rules]];
(* Extend base, if necessary, so that no strong generator fixes all base points *)
If[j===k+1,
AppendTo[newB,Min[Complement[Range[len],StablePoints[gbar,len],newB,{0}]]];
If[verb,Print["Appended B[[",k+1,"]]=",Last[newB]]]
];
(* Ensure we still have a base and SGS for H^(i+1) *)
For[level=j,level≥i+1,--level,{newB,newS}=SchreierSimsStep[newB,newS,len,level,GenSet[gbar],options]];
If[verb,Print["Finished check of H(",i+1,") with base ",Drop[newB,i]," and SGS ",Stabilizer[Take[newB,i],newS]/.rules]];
]

]

];{newB,newS}];

Old computation of j in Butler3:

j=Min[Complement[Range[k],Flatten[Position[newB,#]&/@StablePoints[gbar,len]]]/.{}->{k+1}]

This is the main driver for the Schreier-Sims algorithm. Note the ugly use of the global counter num:

TODO: The resulting base contains a zero if using MathLink->False or Images notation (perhaps other notations too). Detected by Alfonso, March 2011.

(* Main driver *)
Options[SchreierSims]={MathLink:>$xpermQ,UseRules:>{},Method->"Butler3",xPermVerbose->False};
SchreierSims[n_Integer,other__]:=SchreierSims[Range[n],other];
SchreierSims[B_List,Symmetric[inds_List,not_:xAct`xPerm`Cycles],options___]:=StrongGenSet[Sort@inds,TranslatePerm[GenSet@@(xAct`xPerm`Cycles/@Partition[Sort@inds,2,1]),not]];
SchreierSims[B_List,Antisymmetric[inds_List,not_:xAct`xPerm`Cycles],options___]:=StrongGenSet[Sort@inds,TranslatePerm[GenSet@@(-xAct`xPerm`Cycles[#]&/@Partition[Sort@inds,2,1]),not]];
SchreierSims[_,GenSet[],len_Integer,options___]:=StrongGenSet[{},GenSet[]];
SchreierSims[B_List,GS_GenSet,options___?OptionQ]:=SchreierSims[B,GS,PermLength[StrongGenSet[B,GS]],options];
SchreierSims[B_List,GS_GenSet,len_Integer,options:OptionsPattern[]]:=If[OptionValue[MathLink],MathLinkSchreierSims,MathSchreierSims][B,GS,len,options];
SetNumberOfArguments[SchreierSims,{2,Infinity}];
Protect[SchreierSims];
(* Mathematica code *)
MathSchreierSims[B_List,GS_GenSet,len_Integer,options:OptionsPattern[]]:=Module[{k,base=B,PBQ,genset,i,rules,verb},

{rules,verb}=OptionValue[SchreierSims,{options},{UseRules,xPermVerbose}];

num=0;

(* Eliminate identity from GS, if present *)
genset=DeleteCases[GS,ID[First[GS]]];

(* Check initial base *)
base=NonStablePoints[B,genset];
k=Length[base];

(* Recursive call to SchreierSims *)
For[i=k,i≥1,--i,
{base,genset}=Evaluate[SchreierSimsStep[base,genset,len,i,Stabilizer[Take[base,i-1],genset],options]];
];

(* Results *)
If[verb,Print["Checked ",num," generators and obtained group of order ",OrderOfGroup[StrongGenSet[base,genset]]]];
StrongGenSet[DeleteCases[base,0],genset/.rules]

];

SchreierSims[{},GenSet[],4]

RiemannSymmetric:

SchreierSims[{},GenSet[-Perm[{2,1,3,4}],Perm[{3,4,1,2}]],4]

SchreierSims[{},GenSet[Perm[{2,1,3,4,6,5}],Perm[{3,4,1,2,5,6}]],6]

SchreierSims[{},GenSet[-xAct`xPerm`Cycles[{2,1}],xAct`xPerm`Cycles[{1,3},{2,4}]],4]

SchreierSims[{},GenSet[xAct`xPerm`Cycles[{2,1},{5,6}],xAct`xPerm`Cycles[{1,3},{2,4}]],6]

A nontrivial example:

a=xAct`xPerm`Cycles[{1,8,9},{2,11,15},{3,10,12},{4,14,19},{5,16,17},{6,21,20},{7,13,18}];
b=xAct`xPerm`Cycles[{9,18,20},{12,19,17}];
c=xAct`xPerm`Cycles[{10,21,11},{13,16,14}];
d=xAct`xPerm`Cycles[{8,13,21},{10,14,16}];
e=xAct`xPerm`Cycles[{8,16,11},{13,14,21}];
f=xAct`xPerm`Cycles[{2,3,6},{4,7,5}];
g=xAct`xPerm`Cycles[{1,7,6},{3,4,5}];
h=xAct`xPerm`Cycles[{12,20,15},{17,19,18}];

$xPermRules={
a->"a",
b->"b",
c->"c",
d->"d",
e->"e",
f->"f",
g->"g",
h->"h",
Perm[Range[21]]->"id",
TranslatePerm[a,{Perm,21}]->"a",
TranslatePerm[b,{Perm,21}]->"b",
TranslatePerm[c,{Perm,21}]->"c",
TranslatePerm[d,{Perm,21}]->"d",
TranslatePerm[e,{Perm,21}]->"e",
TranslatePerm[f,{Perm,21}]->"f",
TranslatePerm[g,{Perm,21}]->"g",
TranslatePerm[h,{Perm,21}]->"h",
TranslatePerm[a,Rules]->"a",
TranslatePerm[b,Rules]->"b",
TranslatePerm[c,Rules]->"c",
Sort[TranslatePerm[d,Rules]]->"d",
Sort[TranslatePerm[e,Rules]]->"e",
Sort[TranslatePerm[f,Rules]]->"f",
Sort[TranslatePerm[g,Rules]]->"g",
Sort[TranslatePerm[h,Rules]]->"h",
TranslatePerm[a,{Images,21}]->"a",
TranslatePerm[b,{Images,21}]->"b",
TranslatePerm[c,{Images,21}]->"c",
TranslatePerm[d,{Images,21}]->"d",
TranslatePerm[e,{Images,21}]->"e",
TranslatePerm[f,{Images,21}]->"f",
TranslatePerm[g,{Images,21}]->"g",
TranslatePerm[h,{Images,21}]->"h"};

SchreierSims[{},TranslatePerm[GenSet[a,b,c],{Perm,21}],21,xPermVerbose->True,UseRules->$xPermRules]

Note the influence of the notation:

Timing[SchreierSims[{},TranslatePerm[GenSet[a,b,c],{Perm,21}]]]

Timing[SchreierSims[{},GenSet[a,b,c]]]

Timing[SchreierSims[{},TranslatePerm[GenSet[a,b,c],Rules]]]

Timing[SchreierSims[{},TranslatePerm[GenSet[a,b,c],{Images,21}]]]

Another example (Butler p. 133):

a=xAct`xPerm`Cycles[{1,2,4,5,7,3,6}];
b=xAct`xPerm`Cycles[{2,4},{3,5}];
g1=xAct`xPerm`Cycles[{2,6,3,7},{4,5}];
g2=xAct`xPerm`Cycles[{4,7},{5,6}];
g3=xAct`xPerm`Cycles[{4,6},{5,7}];
$xPermRules={a->"a",b->"b",g1->"g1",g2->"g2",g3->"g3"};

SchreierSims[{1,2},GenSet[a,b],UseRules->$xPermRules]

The solution given by Butler is slightly different, because he is using algorithm "Butler2", which we have not encoded. We simply check that his permutation g3 is in our group:

SchreierSims[{1,2},GenSet[a,b]]

PermMemberQ[xAct`xPerm`Cycles[{4,5},{6,7}],%]

4.8. Some important SGS. Symmetric, Antisymmetric, PairSymmetric

Strong generating sets for important groups:

These two are now (from xPerm 1.2.0) inert heads, as Mathematica v9 contains them as well in inert form:

Unprotect[Symmetric];
Symmetric[inds_List,not_]:=StrongGenSet[Sort@inds,TranslatePerm[GenSet@@(Cycles/@Partition[Sort@inds,2,1]),not]];
Symmetric[{},_]:=StrongGenSet[{},GenSet[]];
Protect[Symmetric];

test[SchreierSims,
SchreierSims[{},Symmetric[{5,3,6}]],
StrongGenSet[{3,5,6},GenSet[Cycles[{3,5}],Cycles[{5,6}]]]
]

test[Symmetric,
Symmetric[{5,3,6},Cycles],
StrongGenSet[{3,5,6},GenSet[Cycles[{3,5}],Cycles[{5,6}]]]
]

Unprotect[Antisymmetric];
Antisymmetric[inds_List,not_]:=StrongGenSet[Sort@inds,TranslatePerm[GenSet@@(-Cycles[#]&/@Partition[Sort@inds,2,1]),not]];
Antisymmetric[{},_]:=StrongGenSet[{},GenSet[-ID]];
Protect[Antisymmetric];

test[SchreierSims,
SchreierSims[{},Antisymmetric[{5,3,6}]],
StrongGenSet[{3,5,6},GenSet[-Cycles[{3,5}],-Cycles[{5,6}]]]
]

test[Antisymmetric,
Antisymmetric[{5,3,6},Cycles],
StrongGenSet[{3,5,6},GenSet[-Cycles[{3,5}],-Cycles[{5,6}]]]
]

Here we are including the four points in the base, instead of only the two needed, just for aesthetics of the results when sorting indices:

RiemannSymmetric[inds:{i1_,i2_,i3_,i4_},not_:Cycles]:=StrongGenSet[Sort@inds,TranslatePerm[GenSet[Cycles[{i1,i3},{i2,i4}],-Cycles[{i1,i2}],-Cycles[{i3,i4}]],not]]
RiemannSymmetry=RiemannSymmetric;

addperm[sign:(1|-1)][{p1a_,p1b_},{p2a_,p2b_}]:=GenSet[sign Cycles[{p1a,p2a},{p1b,p2b}]]
addperm[sign:(1|-1)][pair_List,pairs__List]:=GenSet[sign^Length[{pairs}] Cycles@@Transpose@RotateRight[{pair,pairs}],addperm[sign][pairs]]
addperm[_][__]:=GenSet[]

Again, now we add all points in the base:

PairSymmetric[pairs:{{_,_}...},sym1_,sym2_,not_:Cycles]:=StrongGenSet[Sort[Flatten@pairs],TranslatePerm[GenSet[addperm[sym1]@@pairs,Switch[sym2,1,GenSet@@Map[Cycles,pairs],-1,GenSet@@Map[-Cycles[#]&,pairs],_,GenSet[]]],not]];

PairSymmetric[{{1,2},{3,4},{5,6},{7,8},{9,10},{11,12},{13,14}},1,-1]

PairSymmetric[{{1,2},{3,4},{5,6},{7,8},{9,10},{11,12},{13,14}},0,-1]

Note there is no difference in the base between these two results:

PairSymmetric[{{1,2},{3,4}},1,-1]

RiemannSymmetric[{1,2,3,4}]

More examples:

S2:

SGS=SchreierSims[{1},GenSet[Cycles[{1,2}]]]

OrderOfGroup[SGS]

StabilizerChain[SGS]//ColumnForm

S3:

SGS=SchreierSims[{1,2},GenSet[Cycles[{1,2}],Cycles[{2,3}]]]

OrderOfGroup[SGS]

StabilizerChain[SGS]//ColumnForm

S4:

SGS=SchreierSims[{1,2,3},GenSet[Cycles[{1,2}],Cycles[{2,3}],Cycles[{3,4}]]]

OrderOfGroup[SGS]

StabilizerChain[SGS]//ColumnForm

4.9. Stabilization II

4.9.1. Stabilizer

Given a group, this function gives the stabilizer subgroup of a list of points. It first constructs a new SGS with those points at the beginning of the base. Then it is trivial to compute the stabilizer:

Stabilizer[points_List,SGS_]:=Stabilizer[points,SGS,PermLength[SGS]];
Stabilizer[points_List,SGS:StrongGenSet[base_,GS_],len_Integer]:=StrongGenSet[DeleteCases[#1,Alternatives@@points],Stabilizer[points,#2]]&@@BaseChange[SGS,points~Join~DeleteCases[base,Alternatives@@points],len];
SetNumberOfArguments[Stabilizer,{2,3}];
Protect[Stabilizer];

Compare:

SGS=StrongGenSet[{1,2,3},GenSet[Cycles[{1,2}],Cycles[{2,3}]]];

SGS==SchreierSims[SGS[[1]],SGS[[2]],3]

Stabilizer[{2},SGS,3]

Stabilizer[{2},SGS[[2]]]

SGS=.

4.9.2. SetStabilizer

Set-stabilization with strong generating sets is a simple application of the backtrack algorithm. See the Case Study 4 in page 113 of Butler's book.

SetStabilizerProperty[points_List,perm_]:=Complement[OnPoints[points,perm],points]==={};

Options[SetStabilizer]={MathLink:>$xpermQ,xPermVerbose->False,OrderedBase->False};

SetStabilizer[points_,sgs:StrongGenSet[base_,GS_],options:OptionsPattern[]]:=If[OptionValue[MathLink],
MathLinkSetStabilizer[points,Max[points,base,PermLength[GS]],sgs,options],Search[StrongGenSet[base,GS],SetStabilizerProperty[points,#]&,1,StrongGenSet[base,GenSet[]],options]];
SetNumberOfArguments[SetStabilizer,{2,Infinity}];
Protect[SetStabilizer];

Example:

SGS=StrongGenSet[{1,2,3},GenSet[Perm[{2,1,3,4}],Perm[{1,3,2,4}],Perm[{1,2,4,3}]]]

OrderOfGroup[%]

The stabilizer of a set of points is always a subgroup of the set-stabilizer of those points:

Stabilizer[{1,2},SGS]

Dimino[Last[%]]

SetStabilizer[{1,2},SGS]

Dimino[Last[%]]

Another example:

Stabilizer[{1,2,3},SGS]

Dimino[Last[%]]

SetStabilizer[{1,2,3},SGS]

Dimino[Last[%]]

5. Canonicalization

(************************ 5. Canonicalization ***********************)

{TimeUsed[],SessionTime[],TimeUsed[]/SessionTime[]}

5.0. Comments

If[$ReadingVerbose,Print["Reading section 5: Canonicalization"],Null,Null]

The process of canonicalization is based on the use of two groups: the group S of slot-symmetries and the group D of dummy-symmetries. It is natural to send the information concerning S using slots, and the information concerning D using names (that is, positions in the canonical list) of indices, and we shall follow this convention here. However, the papers by Portugal et al. always use slots, having to use a permutation h which keeps track of the shifts of dummies. On the other hand each permutation can be given as it is or through its inverse. Unfortunately we use in xTensor precisely the opposite notation to that of Portugal et al, but still we want to use their algorithms as they are. We solve this conflict by having InversePerm at the input and at the output in the only CanonicalPerm call in xTensor. That means that in this package we shall use Portugal's notation everywhere: permutations map slots to indices: the image of a slot under a permutation is the index at that slot. In other words, sorting according to a basis of S means deciding the priorities of the different slots.

5.1. RightCosetRepresentative

5.1.0. Comments

From the point of view of tensor indices, we must distinguish two different numbered objects: slots and indices. We start from a given index configuration C={a, c, b, ...} which can be obtained from a standard configuration C0={a, b, c, ...} as C.g=C9 with a permutation g. On the other hand we have a symmetry group S of permutations of slots in the tensor, independently of which indices are in those slots. In this sense, two configurations C1 and C2 are equivalent if they can be related by C1=C2.s with s in S. If C1.g1=C0 and C2.g2=C0 then permutations g1 and g2 are equivalent if g2=s.g1. We see that S acts on the left on the permutations (we are working with right-cosets). In other words, s reshuffle the slots and then g maps slots to indices. In the RightCosetRepresentative we always use slots only.

As we said, this is the choice in the Portugal's papers: For instance, in page 5 of math-ph/0107031 he says that permutation (1,3,2,4) in cyclic notation converts T^{abcd} into T^{cdba}. To stress this point, starting version 0.6 (July 2007) of xPerm we have changed the name CosetRepresentative to RightCosetRepresentative.

5.1.1. Code

September 2004. Important change: the subroutine is changed to imitate Canon's behaviour. This is not better or faster, I just want to facilitate comparison of both codes. There are two modifications:
	1) new input: free. List of free slots. Only those slots need to be canonicalized. In our previous code all slots were canonicalized.
	2) new output: apart from the canonical representative, we output the new list of free slots and the new SGS for the group S.

February 2008. The If[Delta1=={},Continue[]] statement requires some additional comments. When this is executed we are leaving a base point without stabilization, and hence we will later stabilize base points which are not first in the base. This is not dangerous because of the following reason: if Delta1 is empty that means that the bi point under analysis belongs to an orbit Delta with no intersection with the free slots. That means that the slots in that orbit are decoupled from the eventual moves of other slots to canonicalize the given permutation with respect to free indices. When we find an orbit in which our free slots participate then that orbit is arranged properly until we drop all free slots from it, and so we go back to the previous case.

?RightCosetRepresentative

Options[xAct`xPerm`RightCosetRepresentative]={
xPermVerbose->False,
BaseChangeCheck->False,
MathLink:>$xpermQ
};
xAct`xPerm`RightCosetRepresentative[p_,len_,sym:(_Symmetric|_Antisymmetric),rest___]:=xAct`xPerm`RightCosetRepresentative[p,len,SchreierSims[{},sym],rest];
xAct`xPerm`RightCosetRepresentative[p_?PermQ,len_Integer,StrongGenSet[base_List,GS_GenSet],free_List,options:OptionsPattern[]]:=Module[{perm=p,newGS=GS,newbase=base,newbase2,newGS2,newfree=free,Delta,Delta1,Deltap,k,pk,pp,om,bi,n=0,verb,bcc},
{verb,bcc}=OptionValue[{xPermVerbose,BaseChangeCheck}];
If[verb,Print["RIGHT-COSET-REPRESENTATIVE ALGORITHM for ",p]];
If[verb,Print["which corresponds to the index list: ",TranslatePerm[InversePerm[p],{Perm,len}]]];
If[verb,Print["base: ",base]];
(* Loop over elements of base *)
For[i=1,i≤Length[base]&&n<Length[free],++i,
bi=base[[i]];
If[verb,Print["****** Analysing element i=",i," of base: slot ",bi," ******"]];
Delta=SchreierOrbit[bi,StrongGenSet[newbase,newGS],len];
If[verb,Print["Symmetry orbit Delta of slots: ",First[Delta]]];
If[verb,Print["Free slots: ",newfree]];
Delta1=Cases[First[Delta],Alternatives@@newfree];
If[verb,Print["Free slots that can go to that slot: ",Delta1]];
If[Delta1=={},Continue[]];
(* Tensor notation: perm acts on slots giving indices *)
Deltap=OnPoints[Delta1,perm];
If[verb,Print["At those slots we respectively find indices Deltap: ",Deltap]];
k=MinB[Deltap,base];
pk=Position[Deltap,k][[1,1]];
If[verb,Print["The least index is ",k,", found at position pk: ",pk," of Deltap"]];
pp=Delta1[[pk]];
If[verb,Print["That index is found in tensor at slot pp: ",pp]];
om=TraceSchreier[pp,Delta];
If[verb,Print["We can move slot ",pp," to slot ",bi," using permutation om: ",om, " in S"]];
perm=PermProduct[om,perm];
If[verb,Print["New indices list: ",TranslatePerm[InversePerm[perm],{Perm,len}]]];
(* New position of free indices *)
newfree=OnPoints[newfree,InversePerm[om]];
If[verb,Print["Computing stabilizer in S of slot ",bi]];
(* Ensure bi is the first point of the base *)
If[verb,Print["SGS before change: ",StrongGenSet[newbase,newGS]]];
{newbase2,newGS2}=List@@BaseChange[StrongGenSet[newbase,newGS],{bi},len];
If[verb,Print["SGS after change: ",StrongGenSet[newbase2,newGS2]]];
newbase2=Drop[newbase2,1];
newGS2=Stabilizer[{bi},newGS2];
If[verb,Print["SGS after stabilization: ",StrongGenSet[newbase2,newGS2]]];
(* Check necessity of the base change *)
If[bcc&&newbase2=!=DeleteCases[newbase,bi],
Print["Checking ",newbase," vs. ", newbase2," with point ",bi];
Block[{$RecursionLimit=Infinity},If[OrderOfGroup[StrongGenSet[newbase2,newGS2]]=!=OrderOfGroup[StrongGenSet[DeleteCases[newbase,bi],Stabilizer[{bi},newGS]]],Print["Unavoidable base change. Contact JMM."]]]];
newbase=newbase2;newGS=newGS2;
If[verb,Print["newbase after change: ",newbase]];
++n;
];
(* We return the canonical permutation and the new SGS *)
{perm,StrongGenSet[newbase,newGS],newfree}
];
(* Default definition *)
xAct`xPerm`RightCosetRepresentative[p_?PermQ,len_Integer,SGS_StrongGenSet,options___]:=xAct`xPerm`RightCosetRepresentative[p,len,SGS,Range[len],options];
SetNumberOfArguments[xAct`xPerm`RightCosetRepresentative,{3,Infinity}];
Protect[xAct`xPerm`RightCosetRepresentative];

5.1.2. Simple examples

Renato's example with the symmetries of a Riemann tensor. Note that he uses base {1,3}: the computation is faster (just two iterations of the loop), but we get an ugly result

reducedSGS=SchreierSims[{},GenSet[-Perm[{2,1,3,4}],Perm[{3,4,1,2}]],4]

SGS=SchreierSims[{1,2,3},GenSet[-Perm[{2,1,3,4}],Perm[{3,4,1,2}]],4]

This permutation represents R[b,c,a,d]:

perm=Perm[{3,1,2,4}];

xAct`xPerm`RightCosetRepresentative[perm,4,reducedSGS,xPermVerbose->True]

%==={-Perm[{1,4,3,2}],StrongGenSet[{},GenSet[]],{4,3,1,2}}

This result is equivalent to -R[a,d,c,b]. Note how the base {1,3} has placed the indices a, c in their canonical positions.

It is better if we fill in the gaps in the base. The computation is slower, but the result is nicer:

xAct`xPerm`RightCosetRepresentative[perm,4,SGS,xPermVerbose->True]

%==={Perm[{1,3,4,2}],StrongGenSet[{},GenSet[]],{3,4,1,2}}

This is equivalent to R[a,d,b,c].

Another permutation in the same right coset (R[d,a,c,b]):

test[xAct`xPerm`RightCosetRepresentative,
xAct`xPerm`RightCosetRepresentative[Perm[{2,4,3,1}],4,SGS],
{Perm[{1,3,4,2}],StrongGenSet[{},GenSet[]],{2,1,4,3}}
]

In other notations:

SGS=TranslatePerm[SGS,xAct`xPerm`Cycles]

test[xAct`xPerm`RightCosetRepresentative,
xAct`xPerm`RightCosetRepresentative[TranslatePerm[perm,xAct`xPerm`Cycles],4,SGS],
{xAct`xPerm`Cycles[{2,4,3}],StrongGenSet[{},GenSet[]],{3,4,1,2}}
]

SGS=TranslatePerm[SGS,Rules]

test[xAct`xPerm`RightCosetRepresentative,
xAct`xPerm`RightCosetRepresentative[TranslatePerm[perm,Rules],4,SGS],
{Rules[2->4,3->2,4->3],StrongGenSet[{},GenSet[]],{3,4,1,2}}
]

SGS=TranslatePerm[SGS,{Images,4}]

perm=TranslatePerm[perm,{Images,4}]

test[xAct`xPerm`RightCosetRepresentative,
xAct`xPerm`RightCosetRepresentative[perm,4,SGS],
{Images[{1,4,2,3}],StrongGenSet[{},GenSet[]],{3,4,1,2}}
]

Let us analyze the situation more closely. This is the right coset S.g, in no particular order:

coset=PermProduct[Dimino[SGS[[2]]],perm]

Using the base {1,3} we see that the least permutation is -Perm[{1,4,3,2}]:

PermSort[coset,{1,3}]

But using the base {1,2,3,4} we see that the least permutation is Perm[{1,4,2,3}]:

PermSort[coset]

We can sort the left coset of inverses and take the inverse again:

coset=PermProduct[InversePerm[perm],Dimino[SGS[[2]]]]

Using the base {1,3} we see that the least permutation is -Perm[{1,4,3,2}], with same inverse

InversePerm/@PermSort[coset,{1,3}]

But using the base {1,2,3,4} we see that the least permutation is Perm[{1,3,4,2}], with inverse Perm[{1,4,2,3}]

InversePerm/@PermSort[coset]

In this case we get the same results using canonicalization of left-cosets and right-cosets. This is not necessarily the case. The problem is that, given a base, the ordering of a set of permutations is not kept if we sort the set of inverses of permutations. For example, take these two permutations which are inverses of each other:

{Images[{3,1,2}],Images[{2,3,1}]}

PermSort[%]

InversePerm/@PermSort[InversePerm/@%%]

In our problem the two alternatives have very different meanings: Sorting left-cosets inv(g).S would mean that we look for images of a base of indices; sorting right-cosets S.g would mean we look for images of a set of slots. The second option is more natural because we have a strong generating set for S, which gives a natural base of slots. There is no natural base of indices apart from the trivial one {1, 2, 3, ...}.

SchreierSims[{},Symmetric[{1,2,3,4,5,6}]]

StrongGenSet[{4,3,2,6,5,1},GenSet[xAct`xPerm`Cycles[{1,2}],xAct`xPerm`Cycles[{2,3}],xAct`xPerm`Cycles[{3,4}],xAct`xPerm`Cycles[{4,5}],xAct`xPerm`Cycles[{5,6}]]]

?SchreierSims

SchreierSims[{4,3,2,6,5,1},GenSet[xAct`xPerm`Cycles[{1,2}],xAct`xPerm`Cycles[{2,3}],xAct`xPerm`Cycles[{3,4}],xAct`xPerm`Cycles[{4,5}],xAct`xPerm`Cycles[{5,6}]],6]

OrderOfGroup[%]

Recall that the list {2,5} gives the initial slots of the free indices, and not the names of the indices (which would be more natural).

xAct`xPerm`RightCosetRepresentative[xAct`xPerm`Cycles[{1,4,6,3,5}],6,Symmetric[{6,5,4,3,1}],{2,5},xPermVerbose->True]

%==={xAct`xPerm`Cycles[{3,4,5,6}],StrongGenSet[{3,4,5,6},GenSet[xAct`xPerm`Cycles[{3,4}],xAct`xPerm`Cycles[{4,5}],xAct`xPerm`Cycles[{5,6}]]],{2,1}}

5.1.3. Nontrivial example. This takes several minutes to evaluate

From the permutations

a=TranslatePerm[xAct`xPerm`Cycles[{1,8,9},{2,11,15},{3,10,12},{4,14,19},{5,16,17},{6,21,20},{7,13,18}],{Images,21}]
b=TranslatePerm[xAct`xPerm`Cycles[{9,18,20},{12,19,17}],{Images,21}]
c=TranslatePerm[xAct`xPerm`Cycles[{10,21,11},{13,16,14}],{Images,21}]

we construct the group

Timing[Length[group=Dimino[GenSet[a,b,c]]]]

and its subgroup

subgroup=Dimino[GenSet[a,b]];
Length[subgroup]

Therefore we have this number of cosets

ncosets=27783/81

SGS=SchreierSims[{},GenSet[a,b],21]

Representatives (perform BaseChange checks):

SetOptions[xAct`xPerm`RightCosetRepresentative,BaseChangeCheck->True];

Timing[result=First[xAct`xPerm`RightCosetRepresentative[#,21,SGS]]&/@group;]

reps=Union[result];
Length[reps]===ncosets

We construct the right cosets of those representatives:

Timing[cosetsR=Table[PermProduct[subgroup,reps[[i]]],{i,ncosets}];]

We check that those cosets form the whole group:

Group@@Sort[Apply[Sequence,cosetsR,1]]===Sort[group]

We also check that we have correctly found all members of all cosets:

Timing[And@@Table[Sort[group[[Flatten@Position[result,reps[[i]]]]]]===Group@@Sort[cosetsR[[i]]],{i,ncosets}]]

Now we check that the representative of each coset is correctly chosen:

Timing[And@@Table[First@PermSort[cosetsR[[i]],{9,8,1}]===reps[[i]],{i,ncosets}]]

Note how different choices give different ways to sort the permutations:

reps[[10]]

coset=List@@cosetsR[[10]];

Note in particular that sorting inverses does not give the same "smallest permutation" than sorting the original set:

PermSort[coset][[Range[5]]]

PermSort[coset,{9,8,1}][[Range[5]]]

InversePerm/@PermSort[InversePerm/@coset,{9,8,1}][[Range[5]]]

5.1.4. Example: S4 in S5

We construct the symmetric group of 5-perms:

s5=Dimino[GenSet[Images[{2,1,3,4,5}],Images[{5,1,2,3,4}]]];
Length[s5]

It is not sorted, according to PermSort:

PermSort[s5]===s5

A subgroup of permutations is s4:

s4=Group@@Stabilizer[{5},GenSet@@s5]

described by

SGSs4=SchreierSims[{},GenSet[Images[{2,1,3,4,5}],Images[{4,1,2,3,5}]],5]

We have 5 coset representatives:

result=First[xAct`xPerm`RightCosetRepresentative[#,5,SGSs4]]&/@List@@s5;

reps=Union[result]

The representatives coincide with the smallest permutations of the right cosets:

coset1=PermSort[PermProduct[s4,reps[[1]]]]

coset1[[1]]==reps[[1]]

coset2=PermSort[PermProduct[s4,reps[[2]]]]

coset2[[1]]==reps[[2]]

coset3=PermSort[PermProduct[s4,reps[[3]]]]

coset3[[1]]==reps[[3]]

coset4=PermSort[PermProduct[s4,reps[[4]]]]

coset4[[1]]==reps[[4]]

coset5=PermSort[PermProduct[s4,reps[[5]]]]

coset5[[1]]==reps[[5]]

Group@@Union[Coset@@coset1,coset2,coset3,coset4,coset5]===Sort[s5]

5.1.5. Example: S3 in S5

We construct the symmetric group of 5-perms:

s5=Dimino[GenSet[Images[{2,1,3,4,5}],Images[{5,1,2,3,4}]]];
Length[s5]

It is not sorted, according to PermSort:

PermSort[s5]===s5

A subgroup of permutations is s3:

s3=Group@@Stabilizer[{4,5},GenSet@@s5]

described by

SGSs3=SchreierSims[{},GenSet[Images[{2,1,3,4,5}],Images[{3,1,2,4,5}]],5]

We have 20 coset representatives:

result=First[xAct`xPerm`RightCosetRepresentative[#,5,SGSs3]]&/@List@@s5;

reps=Union[result]

(cosets=Table[PermSort[PermProduct[s3,reps[[i]]]],{i,1,20}])//ColumnForm

And@@Table[cosets[[i,1]]==reps[[i]],{i,1,20}]

Group@@Union[Apply[Sequence,cosets,1]]===Sort[s5]

5.1.6. Example: S3 in S6

We construct the symmetric group of 6-perms:

s6=Dimino[GenSet[Images[{2,1,3,4,5,6}],Images[{6,1,2,3,4,5}]]];
Length[s6]

It is not sorted, according to PermSort:

PermSort[s6]===s6

A subgroup of permutations is s3:

s3=Group@@Stabilizer[{4,5,6},GenSet@@s6]

described by

SGSs3=SchreierSims[{},GenSet[Images[{2,1,3,4,5,6}],Images[{3,1,2,4,5,6}]],6]

We have 20 coset representatives:

result=First[xAct`xPerm`RightCosetRepresentative[#,6,SGSs3]]&/@List@@s6;

reps=Union[result];

cosets=Table[PermSort[PermProduct[s3,reps[[i]]]],{i,1,120}];

And@@(Table[cosets[[i,1]]==reps[[i]],{i,1,120}])

Group@@Union[Apply[Sequence,cosets,1]]===Sort[s6]

5.2. DoubleCosetRepresentative

5.2.1. DummySet and RepeatedSet

We first define the concept of a dummyset. We need it because we can have a tensor with dummies coming from different manifolds. A DummySet is a set of pairs of dummies all belonging to the same manifold. In principle there is no reason to work simultaneously with two disjoint dummy sets on the same manifold. Therefore we also require that a DummySet is the set of all pairs of dummies belonging to the same manifold. In addition, we supply the symmetry of the metric of that manifold, with this convention: if the metric is symmetric then metricQ=1; if the metric is antisymmetric (spinor calculus) then metricQ=-1; if there is no metric in the manifold then metricQ=0. Note that in the list of pairs of dummies we are giving their names, that is, we are supplying their positions in the canonical configuration. This is because D is a group acting on indices and not on slots.

The DoubleCosetAlgorithm requires changing the SGS of D by BaseChange and stabilization. However, it is easier to change the dummysets and reconstruct the associated SGS because it always has the same simple structure. This is done with SGSOfDummySet:

?SGSOfDummySet

SGSOfDummySet[DummySet[_,pairs_List,metricQ_Integer]]:=StrongGenSet[First/@pairs,
GenSet@@Join[Switch[metricQ,
1,Cycles/@pairs,
-1,Minus/@Cycles/@pairs,
0,{},
_,Throw[Print["Invalid value for metricQ in SGSOfDummySet."]]
],
Flatten[Cycles@@#&/@Transpose/@Partition[pairs,2,1]]]];
SetNumberOfArguments[SGSOfDummySet,1];
Protect[DummySet];

There are two operations defined on a dummyset: moving and removing indices. Removing a pair is straightforward:

RemovePairOf[i_Integer,DummySet[m_,pairs_List,metricQ_Integer]]:=DummySet[m,DeleteCases[pairs,{___,i,___}],metricQ];

However, moving an index to the first position requires the convention of having all first indices with the same character and all second indices with the same character. This is because the exchange of dummies must be done among up-indices and among down-indices separately. It might seem that we need to worry about a sign in the spinor case in the second line, but this is not true: the actual order of the pairs is not relevant and that sign is stored in the permutations d1.

MovePairOf[i_Integer,DummySet[m_,{p1___,{i_,j_},p2___},metricQ_Integer]]:=DummySet[m,{{i,j},p1,p2},metricQ];
MovePairOf[i_Integer,DummySet[m_,{p1___,{j_,i_},p2___},metricQ_Integer]]:=DummySet[m,Reverse/@{{j,i},p1,p2},metricQ];
MovePairOf[i_Integer,DS_DummySet]:=DS;

Examples:

SGSOfDummySet[DummySet[M2,{{1,2},{3,4}},0]]

SGSOfDummySet[DummySet[M2,{{1,2}},1]]

SGSOfDummySet[DummySet[M2,{{1,2},{3,4}},1]]

SGSOfDummySet[DummySet[M2,{{1,2},{3,4},{5,6}},1]]

SGSOfDummySet[DummySet[M2,{{1,2},{3,4}},-1]]

RemovePairOf[4,DummySet[M2,{{1,2},{3,4},{5,6}},1]]

MovePairOf[3,DummySet[M2,{{1,2},{3,4},{5,6}},1]]

MovePairOf[4,DummySet[M2,{{1,2},{3,4},{5,6}},1]]

MovePairOf[7,DummySet[M2,{{1,2},{3,4},{5,6}},1]]

Then we imitate a dummyset using the (simpler) structure of a repeatedset. It only contains a list of points, but that list can have arbitrary length (if it has length 0 then it is dropped from the list of repeatedsets). We allow cases with length 1 because this will be used in xCoba for nicer canonicalization of c-indices.

RepeatedSet[{}]:=Sequence[];
SGSOfDummySet[RepeatedSet[list_List]]:=SchreierSims[{},Symmetric[list]];
SGSOfDummySet[list_List]:=Apply[JoinSGS,SGSOfDummySet/@list];
MovePairOf[i_Integer,RepeatedSet[{p1___,i_,p2___}]]:=RepeatedSet[{i,p1,p2}];
MovePairOf[i_Integer,repe_RepeatedSet]:=repe;
RemovePairOf[i_Integer,RepeatedSet[{p1___,i_,p2___}]]:=RepeatedSet[{p1,p2}];
RemovePairOf[i_Integer,repe_RepeatedSet]:=repe;
Protect[SGSOfDummySet,RepeatedSet];

RepeatedSet[{}]

RepeatedSet[{4}]

SGSOfDummySet[RepeatedSet[{1,3,5}]]

MovePairOf[3,RepeatedSet[{1,2,3,4,5}]]

MovePairOf[6,RepeatedSet[{1,2,3,4,5}]]

RemovePairOf[3,RepeatedSet[{1,2,3,4,5}]]

RemovePairOf[6,RepeatedSet[{1,2,3,4,5}]]

5.2.2. JoinSGS

Joining of disjoint Strong Generating Sets (we assume that the points are already correctly numbered). This can be understood as the construction of the SGS for the direct product of several groups.

?JoinSGS

JoinSGS[prev___,sym:(_Symmetric|_Antisymmetric),after___]:=JoinSGS[prev,SchreierSims[{},sym],after];
JoinSGS[SGSs__StrongGenSet]:=StrongGenSet[Join@@#1,Union@@#2]&@@Transpose[Apply[List,{SGSs},1]];
JoinSGS[]:=StrongGenSet[{},GenSet[]];

JoinSGS[StrongGenSet[{1,3},GenSet[A,B,C]],StrongGenSet[{5,6},GenSet[D,F]]]

5.2.3. ExtendBase

The following function extends the original base to include all dummyslots, for sorting purposes, but does not change the relative ordering of the initial points in the base. This algorithm can be used in a general case, but it was designed for use in DoubleCosetRepresentative calling from CanonicalPerm, such that all points in base are dummies. We check that and print a message when this is not true.

takenext[{a_,other___},{found___},list_List]:=With[{sublist=Complement[Select[list,#<a&],{other}]},takenext[{other},Flatten[{found,sublist,a}],Complement[list,sublist,{a}]]];
takenext[{},found_,list_List]:=Join[found,list];
ExtendBase[base_,dummyslots_,ob_]:=Module[{nondummybase=Complement[base,dummyslots],otherdummyslots=Complement[dummyslots,base]},
If[nondummybase=!={},Print["ExtendBase: non-dummy base points: ",nondummybase]];
If[ob,
takenext[base,{},dummyslots],
Join[Cases[base,Alternatives@@dummyslots],otherdummyslots]
]
];

TracePrint[takenext[{1,5,3,7},{},{1,2,3,4,5,6,7}],takenext[_List,_List,_List]]

This is here only for this notebook (the full list of options for CanonicalPerm will be given below, and it includes this in particular).

Options[CanonicalPerm]:={OrderedBase->True};

5.2.4. Main code

Simple function:

orbitof[i_Integer,h_[___,list:{___,i_,___},___]]:=list
orbitof[i_Integer,_[___]]:={}

List of objects in the code:
permutation:	permutation to canonicalize
perm:		permutation with notation ID corrected.
notation:	inner notation to work with
dummyindices: indices (positions in the canonical list) participating in the D group
dummyslots:	list obtained as the union of all positions (that is, slots: positions in the actual configuration) participating in the D group

bS:		base of the SGS of group S. It must contain all positions in dummyslots, and so base is extended if this is not true initially. The OrderedBase option add missing dummyslots in (Sort) order, without moving the original slots to avoid losing the base property.
SGSS:		SGS of group S, initially formed by bS and the given GS. This GS is assumed to have been stabilized on the free (and blocked) points.
bSsort:	rearrangement for sorting. This is in fact a choice of ordering of indices.

SGSD:		SGS of group D, initially formed by joining all SGS constructed with SGSOfDummySet.

Now this function finds the double-coset representative. The degree deg that is given as an argument is the the degree of the permutation g, and does not have to coincide with the number of dummies.

?DoubleCosetRepresentative

Options[DoubleCosetRepresentative]:={
xPermVerbose->False,
MathLink:>$xpermQ
};
(* 
permutation: degree deg;
base={b1, ..., bk}: base of group S ;
GS: strong generating set of S relative to base ;
dummysets : list of sets (head DummySet or RepeatedSet) of pairs of names
(positions in the canonical configuration) of dummies, or a single list
with repeated indices. Each set is associated with a manifold ;
*)
DoubleCosetRepresentative[permutation_?PermQ,len_Integer,StrongGenSet[base_List,GS_GenSet],dummysets:{(_DummySet|_RepeatedSet)...},options:OptionsPattern[]]:=Module[{n,perm,notation,F1,F2,bS,SGSS,SGSD,i,TAB,ALPHA={{}},nuS,Deltab,DeltaD,IMAGES,p,nuD,Deltap,s,d,NEXT,j,jj,s1,d1,L1,tmp,alphaindices,KS=GS,result,dummyindices,newdummysets=dummysets,dummyslots,bi,pi,ob=OptionValue[CanonicalPerm,{options},OrderedBase],verb=OptionValue[xPermVerbose],bSsort},

(* 1. The routine works in the notation given by permutation. If it is ID, change to Cycles[] *)
perm=permutation/.ID->xAct`xPerm`Cycles[];
(* Note that perm is always kept fixed *)
notation=NotationOfPerm[perm];
If[verb,Print["DOUBLE-COSET-REPRESENTATIVE ALGORITHM for ",permutation]];

(* 2. Dummies and repeated indices go into the D group *)
If[verb,Print["index-dummysets: ",dummysets]];
dummyindices=Flatten[dummysets/.{DummySet[_,dums_,_]->dums,RepeatedSet[list_]->list}];
If[verb,Print["dummyindices: ",dummyindices]];
dummyslots=OnPoints[dummyindices,InversePerm[perm]];
If[verb,Print["dummyslots: ",dummyslots]];
newdummysets=dummysets;

(* 3. Extend base. base contains some of the dummy slots, but not all. We must cover all slots of dummies, keeping the order of the points of the original base *)
bS=ExtendBase[base,Sort@dummyslots,ob];
If[verb,Print["Extended ",base," to ",bS]];
SGSS=StrongGenSet[bS,GS];
If[verb,Print["Initial SGSS: ",SGSS]];
(* Initialize base images *)
p=0 bS;

(* 4. Adapt base for sorting of dummy indices. This is a particular choice, for aesthetical reasons *)
bSsort=bS/.Inner[Rule,Sort[bS],Sort[dummyindices],List];
If[verb,Print["base for sorting: ",bSsort]];

(* 5. Initialize TAB *)
TAB[{}]:={ID[perm],ID[perm]};

(* 6. Subroutines. Note that elements of TAB are {s, d} *)
F2[L_,TAB_,g_]:=PermProduct[TAB[L][[1]],g,TAB[L][[2]]];
F1[L_,TAB_,DeltaD_,Deltab_,g_]:=Module[{sgd=F2[L,TAB,g],list},
If[verb,Print["With L=",L," we get sgd: ",sgd]];
list=OnPoints[Deltab,sgd];
If[verb,Print["which maps slots in Deltab to indices list: ",list]];
list=Union[Flatten[Cases[DeltaD,x_/;Intersection[x,list]=!={}]]];
If[verb,Print["whose points belong to orbits ",list]];
list
];

(* 7. Strong Generating Set of group D *)
SGSD=TranslatePerm[JoinSGS@@(SGSOfDummySet/@dummysets),notation];
If[verb,Print["Initial SGSD: ",SGSD]];

(* 8. Main loop on the slots of bS *)
For[i=1,i≤Length[bS],++i,
If[verb,Print["******************* Loop i= ",i," *********************"]];
bi=bS[[i]];
If[verb,Print["Analyzing slot ",bi," of tensor"]];

(* A. Schreier vector of S *)
nuS=SchreierOrbits[SGSS,len];
(* Orbits are ordered according to base of SGSS, which has bi as first element *)
If[verb,Print["nuS: ",nuS," with first element ",bi]];
Deltab=orbitof[bi,Drop[List@@nuS,-2]];
If[verb,Print["Under S, slot ",bi," can go to slots Deltab: ",Deltab]];

(* B. Orbits of D *)
DeltaD=Orbits[SGSD[[2]],len];
If[verb,Print["Orbits of indices under D: DeltaD: ",DeltaD]];

(* C. Images of bi under elements of S.perm.D *)
IMAGES=Union[Flatten[Map[F1[#,TAB,DeltaD,Deltab,perm]&,ALPHA]]];
If[verb,Print["Therefore at slot ",bi," we can have indices IMAGES: ",IMAGES]];
(* The minimal element is taken with respect to bSsort *)
p[[i]]=pi=MinB[IMAGES,bSsort];
If[verb,Print["The least of them is p[[",i,"]]: ",pi]];

(* D. Rearrange SGS of D. Now pi will be the first element of the base *)
newdummysets=MovePairOf[pi,#]&/@newdummysets;
If[verb,Print["Moved pairs ",newdummysets]];
(* Full reconstruction of the SGSD starting from the new dummysets, with reordered indices *)
SGSD=TranslatePerm[JoinSGS@@(SGSOfDummySet/@newdummysets),notation];
If[verb,Print["New SGS of D: ",SGSD]];
(* Schreier vector of D *)
nuD=SchreierOrbits[SGSD,len];
If[verb,Print["with Schreier vector nuD: ",nuD]];
(* Orbit of pi under D *)
Deltap=orbitof[pi,Drop[List@@nuD,-2]];
If[verb,Print["In particular, the orbit of index ",pi," is Deltap: ",Deltap]];

(* E. Calculate ALPHA and TAB. Inner double loop *)
If[verb,Print["Now looking for all permutations sgd that move index ",pi," to slot ",bi]];
alphaindices={};
For[l=1,l≤Length[ALPHA],++l,
If[verb,Print["Loop with l=",l]];
L=ALPHA[[l]];
If[verb,Print["L= ",L]];
s=TAB[L][[1]];
d=TAB[L][[2]];
If[verb,Print["TAB[L] = {s, d} = ",{s,d}]];
If[verb,Print["Calculating NEXT. We need the intersection of sets of slots ", OnPoints[Deltab,s]," and ",OnPoints[Deltap,InversePerm[PermProduct[perm,d]]]]];
NEXT=Intersection[OnPoints[Deltab,s],OnPoints[Deltap,InversePerm[PermProduct[perm,d]]]];
If[verb,Print["Intermediate slots NEXT= ",NEXT]];
For[jj=1,jj≤Length[NEXT],++jj,
j=NEXT[[jj]];
s1=PermProduct[TraceSchreier[OnPoints[j,InversePerm@s],nuS],s];
If[verb,Print["From slot ",bi," to intermediate slot ",j," use s1=", s1]];
d1=PermProduct[d,InversePerm@TraceSchreier[OnPoints[j,PermProduct[perm,d]],nuD]];
If[verb,Print["d1= ",d1]];
L1=Append[L,j];
If[verb,Print["L1= ",L1]];
TAB[L1]={s1,d1};
AppendTo[alphaindices,L1];
If[verb,Print["This gives us the new index configuration: ",TranslatePerm[InversePerm[F2[L1,TAB,perm]],{Perm,len}]]];
(* Checks *)
For[ii=1,ii≤i,++ii,If[OnPoints[bS[[ii]],PermProduct[s1,perm,d1]]=!=p[[ii]],Print["WRONG check of slot ",bS[[ii]]," with point ",p[[ii]]]]
]];
TAB[L]=.;
];
ALPHA=alphaindices;
If[verb,Print["New ALPHA: ",ALPHA]];

(* F. Verify if there are 2 equal permutations of opposite sign in S.perm.D *)
tmp=Map[InversePerm[F2[#,TAB,perm]]&,ALPHA];
If[verb,Print["Checking consistency in set ",tmp]];
If[Intersection[tmp,-tmp]=!={},result=0;Break[]];

(* G. Find the stabilizers S^(i+1) and D^(i+1) *)
SGSS=StrongGenSet[Drop[SGSS[[1]],1],Stabilizer[{bi},SGSS[[2]]]];
If[verb,Print["Removing permutations from SGS of S that move slot ",bi]];
If[verb,Print["New SGS of S: ",SGSS]];
(* Construct new SGSD from the new dummysets, and not by stabilization *)
newdummysets=RemovePairOf[pi,#]&/@newdummysets;
SGSD=TranslatePerm[JoinSGS@@(SGSOfDummySet/@newdummysets),notation];
If[verb,Print["Removing permutations from SGS of D that move index ",pi]];
If[verb,Print["New SGS of D: ",SGSD]];
]; (* End of main loop *);

(* 9. Result *)
Switch[result,
0,0,
_,F2[ALPHA[[1]],TAB,perm]/.If[permutation===ID,Cycles[]->ID,{}]
]

];
SetNumberOfArguments[DoubleCosetRepresentative,{4,Infinity}];
Protect[DoubleCosetRepresentative];

5.2.5. Simple examples

Two-index tensors with both indices contracted:

emptySGS=StrongGenSet[{},GenSet[]];

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{1,2}],2,emptySGS,{DummySet[M,{{1,2}},1]}],
Perm[{1,2}]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Cycles[],2,emptySGS,{DummySet[M,{{1,2}},1]}],
Cycles[]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Rules[],2,emptySGS,{DummySet[M,{{1,2}},1]}],
Rules[]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[ID,2,emptySGS,{DummySet[M,{{1,2}},1]}],
ID
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Images[{1,2}],2,emptySGS,{DummySet[M,{{1,2}},1]}],
Images[{1,2}]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{2,1}],2,emptySGS,{DummySet[M,{{1,2}},1]}],
Perm[{1,2}]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{2,1}],2,emptySGS,{DummySet[M,{{1,2}},1]}],
Perm[{1,2}]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Cycles[{1,2}],2,emptySGS,{DummySet[M,{{1,2}},1]}],
Cycles[]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Rules[1->2,2->1],2,emptySGS,{DummySet[M,{{1,2}},1]}],
Rules[]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Cycles[],4,emptySGS,{DummySet[M,{{3,4}},1]}],
Cycles[]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Cycles[],4,emptySGS,{DummySet[M,{{1,2},{3,4}},1]}],
Cycles[]
]

Now the tensors are symmetric or antisymmetric:

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{1,2}],2,StrongGenSet[{1},GenSet[Perm[{2,1}]]],{DummySet[M,{{1,2}},1]}],
Perm[{1,2}]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{1,2}],2,StrongGenSet[{1},GenSet[-Perm[{2,1}]]],{DummySet[M,{{1,2}},1]}],
0
]

Four-indices tensors:

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{1,2,3,4}],4,emptySGS,{DummySet[M,{{1,2},{3,4}},1]}],
Perm[{1,2,3,4}]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{3,4,1,2}],4,emptySGS,{DummySet[M,{{1,2},{3,4}},1]}],
Perm[{1,2,3,4}]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{1,2,3,4}],4,StrongGenSet[{1},GenSet[-Perm[{2,1,3,4}]]],{DummySet[M,{{1,2},{3,4}},1]}],
0
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{3,4,1,2}],4,StrongGenSet[{3},GenSet[-Perm[{1,2,4,3}]]],{DummySet[M,{{1,2},{3,4}},1]}],
0
]

This is an example of a case in which we have not stabilized the free points. There is a message complaining about point 3 being lost.

Catch@DoubleCosetRepresentative[Perm[{3,4,2,1}],4,StrongGenSet[{3},GenSet[-Perm[{1,2,4,3}]]],{DummySet[M,{{3,4}},1]}]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{3,4,1,2}],4,StrongGenSet[{1},GenSet[-Perm[{2,1,3,4}]]],{DummySet[M,{{3,4}},1]}],
0
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{3,4,1,2}],4,StrongGenSet[{1,2,3},GenSet[Perm[{1,2,4,3}]]],{DummySet[M,{{1,2},{3,4}},1]}],
Perm[{1,2,3,4}]
]

This case is T^BA_BA with symmetry in the third and fourth indices, and we get T^AB_AB

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{2,4,1,3}],4,StrongGenSet[{1,2,3},GenSet[Perm[{1,2,4,3}]]],{DummySet[M,{{1,2},{3,4}},1]}],
Perm[{1,3,2,4}]
]

This case is T^AB_BA with antisymmetry in the third and fourth indices, and we get -T^AB_AB

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Perm[{1,4,2,3}],4,StrongGenSet[{1},GenSet[-Perm[{1,2,4,3}]]],{DummySet[M,{{1,2},{3,4}},1]}],
-Perm[{1,3,2,4}]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Cycles[{2,3}],4,StrongGenSet[{},GenSet[-Cycles[{2,3}]]],{DummySet[M,{{1,2},{3,4}},1]}],
-Cycles[]
]

test[DoubleCosetRepresentative,
DoubleCosetRepresentative[Rules[2->3,3->2],4,StrongGenSet[{1},GenSet[-Rules[2->3,3->2]]],{DummySet[M,{{1,2},{3,4}},1]}],
-Rules[]
]

And@@alltests[DoubleCosetRepresentative]

5.2.6. Riemann tensor

In the past the DoubleCosetRepresentative algorithm was able to canonicalize also free indices. Now it is not, but it is faster because it only worries about dummies.

All configurations:

configurations=Perm@#&/@Permutations[{1,2,3,4}]

SGS=StrongGenSet[{1,3},GenSet[Perm[{3,4,1,2}],-Perm[{2,1,3,4}],-Perm[{1,2,4,3}]]];
result=DoubleCosetRepresentative[#,4,SGS,{DummySet[M,{{1,2},{3,4}},1]}]&/@configurations

Union[result]

which are R^{ab}_{ab} and its negative.

5.3. Both together. CanonicalPerm

5.3.1. Code

INPUT:
	perm: signed permutation in any format
	len: degree of permutation perm
	GS: generating set of a subgroup of S_len (numbered slots)
	frees: free indices (their positions in the canonical list)
	dummysets: dummies and repeated indices (their positions in the canonical list)

length[dummysets_List]:=Plus@@(length/@dummysets);
length[DummySet[_,pairs_List,_]]:=2Length[pairs];
length[RepeatedSet[list_List]]:=Length[list];

?CanonicalPerm

Note the option OrderedBase->True for nicer output. Thomas has found (June 2011) that results can be incorrect with OrderedBase->True if the base does not have all moved points and they are sorted. We assume such a thing now in all of xAct.

Options[CanonicalPerm]={MathLink:>$xpermQ,TimeVerbose->False,xPermVerbose->False,OrderedBase->True};
CanonicalPerm[perm_?PermQ,len_Integer,GS_,frees_List,dummysets:{(_DummySet|_RepeatedSet)...},options:OptionsPattern[]]:=If[OptionValue[MathLink],MathLinkCanonicalPerm,MathCanonicalPerm][perm,len,GS,frees,dummysets,options];
SetNumberOfArguments[CanonicalPerm,{5,Infinity}];
Protect[CanonicalPerm];

MathCanonicalPerm[perm_,len_,GS_GenSet,frees_,dummysets_,options:OptionsPattern[]]:=Module[{notation,SGS,time=TimeUsed[],timeverb,verb,ob},

{timeverb,verb,ob}=OptionValue[CanonicalPerm,{options},{TimeVerbose,xPermVerbose,OrderedBase}];

(* Compute Strong Generating Set *)
notation=NotationOfPerm[perm];
If[verb,Print["Using notation ",notation]];
SGS=Sequence[If[ob,Range[len],NonStablePoints[{},GS]],TranslatePerm[GS,notation]];
If[verb,Print["Before SchreierSims: ",{SGS}]];
SGS=SchreierSims[SGS,len,options];
If[verb,Print["After SchreierSims: ",SGS]];

If[timeverb,Print["SGS for group of order ",OrderOfGroup[SGS]," computed in ",-time+(time=TimeUsed[])," secs."]];

(* {perm, frees, dummysets} come directly here *)
MathCanonicalPerm[perm,len,SGS,frees,dummysets,options]
];

MathCanonicalPerm[perm_,len_,sym:(_Symmetric|_Antisymmetric),rest__]:=MathCanonicalPerm[perm,len,SchreierSims[{},sym],rest];
MathCanonicalPerm[perm_,len_,SGS_StrongGenSet,frees_,dummysets_,options:OptionsPattern[]]:=Module[{notation,freeps,newperm,newSGS,newfreeps,time=TimeUsed[],timeverb,verb,ob},

{timeverb,verb,ob}=OptionValue[CanonicalPerm,{options},{TimeVerbose,xPermVerbose,OrderedBase}];

(* Length of permutations (without considering the sign) and notation *)
If[len=!=length[dummysets]+Length[frees],Throw@Print["Wrong permutation length."]];
notation=NotationOfPerm[perm];

(* Strong generating set *)
newSGS=TranslatePerm[SGS,notation];
If[ob,newSGS=ReplacePart[newSGS,takenext[First@newSGS,{},Range[len]],1]];

(* First, we apply the free-indices algorithm. Compute slots of frees *)
freeps=OnPoints[frees,InversePerm[perm]];
If[verb,Print["Free indices at slots: ",freeps]];
{newperm,newSGS,newfreeps}=RightCosetRepresentative[perm,len,newSGS,freeps,options];
If[timeverb,Print["Free algorithm applied in ",-time+(time=TimeUsed[])," secs."]];
If[verb,Print["Canonical Permutation after RightCosetRepresentative: ",newperm]];

(* If there are dummies apply dummy-indices algorithm*)
If[dummysets=!={},
newperm=DoubleCosetRepresentative[newperm,len,newSGS,dummysets,options];
If[timeverb,Print["Dummy algorithm applied in ",-time+TimeUsed[]," secs."]];
];

(* Final result (could be zero) *)
newperm
];

5.3.2. Examples

Use the Mathematica code only, in this section:

SetOptions[CanonicalPerm,MathLink->False]

SGS=StrongGenSet[{1,3,5,7,9,11,13,15,17,19},GenSet[Cycles[{9,13},{10,14},{11,15},{12,16}],Cycles[{9,13},{10,14},{11,15},{12,16}],-Cycles[{19,20}],Cycles[{17,19},{18,20}],-Cycles[{17,18}],-Cycles[{15,16}],Cycles[{13,15},{14,16}],-Cycles[{13,14}],-Cycles[{11,12}],Cycles[{9,11},{10,12}],-Cycles[{9,10}],-Cycles[{7,8}],Cycles[{5,7},{6,8}],-Cycles[{5,6}],Cycles[{13,17},{14,18},{15,19},{16,20}],Cycles[{1,5},{2,6},{3,7},{4,8}],-Cycles[{1,2}],Cycles[{1,3},{2,4}],-Cycles[{3,4}]]];

test[CanonicalPerm,
CanonicalPerm[Cycles[{2,10,12,17,20,18,11,6,14,19,7,4,13,8,15,9,5,3}],20,SGS,{},{DummySet[M,{{1,2},{3,4},{5,6},{7,8},{9,10},{11,12},{13,14},{15,16},{17,18},{19,20}},1]},TimeVerbose->True],
0
]

test[CanonicalPerm,
CanonicalPerm[Perm[{3,1,2,4}],4,GenSet[-Perm[{2,1,3,4}],Perm[{3,4,1,2}]],{1,2,3,4},{},TimeVerbose->True,OrderedBase->False],
-Perm[{1,4,3,2}]
]

test[CanonicalPerm,
CanonicalPerm[Perm[{3,1,2,4}],4,GenSet[-Perm[{2,1,3,4}],Perm[{3,4,1,2}]],{1,2,3,4},{},TimeVerbose->True],
Perm[{1,3,4,2}]
]

Renato's example in paper II:

TranslatePerm[Cycles[{1,12,11,4,5,6,8,9,10,7,3}],Perm]

sperm=Perm[{3,2,7,11,4,5,10,6,8,9,12,1}];
GS=GenSet[-Perm[{2,1,3,4,5,6,7,8,9,10,11,12}],-Perm[{1,2,4,3,5,6,7,8,9,10,11,12}],-Perm[{1,2,3,4,6,5,7,8,9,10,11,12}],-Perm[{1,2,3,4,5,6,8,7,9,10,11,12}],-Perm[{1,2,3,4,5,6,7,8,10,9,11,12}],-Perm[{1,2,3,4,5,6,7,8,9,10,12,11}],Perm[{3,4,1,2,5,6,7,8,9,10,11,12}],Perm[{1,2,3,4,7,8,5,6,9,10,11,12}],Perm[{1,2,3,4,5,6,7,8,11,12,9,10}],Perm[{1,2,3,4,9,10,11,12,5,6,7,8}]];

Here all indices are free:

test[CanonicalPerm,
CanonicalPerm[sperm,12,GS,Range[1,12],{},TimeVerbose->True],
Perm[{1,3,5,9,2,7,11,8,6,12,10,4}]
]

test[CanonicalPerm,
CanonicalPerm[sperm,12,GS,Range[1,12],{},TimeVerbose->True,OrderedBase->False],
-Perm[{1,3,5,12,2,7,9,8,6,10,11,4}]
]

Here we include the dummy indices:

test[CanonicalPerm,
CanonicalPerm[sperm,12,GS,{1,2},{DummySet[M,{{3,4},{5,6},{7,8},{9,10},{11,12}},1]},TimeVerbose->True],
-Perm[{1,3,2,5,4,9,6,11,7,10,8,12}]
]

and this is the answer given by Renato:

TranslatePerm[-Perm[{1,3,2,5,4,9,6,11,7,10,8,12}],Cycles]

-System`Cycles[{{3,2},{5,4},{7,9,6},{11,8}}]

Without ordering the base we get a different result:

test[CanonicalPerm,
CanonicalPerm[sperm,12,GS,{1,2},{DummySet[M,{{3,4},{5,6},{7,8},{9,10},{11,12}},1]},TimeVerbose->True,OrderedBase->False],
-Perm[{1,3,2,12,5,4,7,10,9,6,11,8}]
]

TranslatePerm[-Perm[{1,3,2,12,5,4,7,10,9,6,11,8}],Cycles]

-System`Cycles[{{3,2},{6,10,8,12,4}}]

test[CanonicalPerm,
CanonicalPerm[TranslatePerm[sperm,Rules],12,TranslatePerm[GS,Rules],{1,2},{DummySet[M,{{3,4},{5,6},{7,8},{9,10},{11,12}},1]},TimeVerbose->True],
-Rules[2->3,3->2,4->5,5->4,6->7,7->9,8->11,9->6,11->8]
]

A trivial case:

test[CanonicalPerm,
CanonicalPerm[Perm[{1,2,3,4}],4,GenSet[],{},{DummySet[M,{{1,2},{3,4}},1]}],
Perm[{1,2,3,4}]
]

More examples:

test[CanonicalPerm,
CanonicalPerm[Perm[{6,5,4,3,2,1}],6,GenSet[-Cycles[{3,5}],-Cycles[{5,6}]],{1,2,3,4,5,6},{},OrderedBase->True],
-Perm[{3,5,4,6,2,1}]
]

Timing[CanonicalPerm[Perm[{6,5,4,3,2,1}],6,GenSet[-Cycles[{3,5}],-Cycles[{5,6}]],{1,2,3,4,5,6},{}]]

Timing[CanonicalPerm[Perm[{6,5,4,3,2,1}],6,Antisymmetric[{5,3,6}],{1,2,3,4,5,6},{}]]

Timing[CanonicalPerm[Perm[{3,2,1,4}],4,RiemannSymmetric[{1,2,3,4}],{1,2,3,4},{}]]

Note the differences in timings with respect to the permuted points and the option OrderedBase:

Timing[CanonicalPerm[Cycles[{1,3}],4,RiemannSymmetric[{1,2,3,4}],Range[4],{}]]

Timing[CanonicalPerm[Cycles[{1,3}],4,RiemannSymmetric[{1,2,3,4}],Range[4],{},OrderedBase->False]]

Timing[CanonicalPerm[Cycles[{1001,1003}],1004,RiemannSymmetric[{1001,1002,1003,1004}],Range[1004],{}]]

Timing[CanonicalPerm[Cycles[{1001,1003}],1004,RiemannSymmetric[{1001,1002,1003,1004}],Range[1004],{},OrderedBase->False]]

The point is that the base in the OrderedBase->True case has 1004 points, but the base in the OrderedBase->False case has only 2 points!

5.4. Transversals

5.4.1. Right coset transversals

Given a SGS of a subgroup of S[n], we want to construct a transversal of representatives of the right cosets. We do it by generation of random permutations.

nosign[-x_]:=x;
nosign[x_]:=x;

?RightTransversal

RightTransversal[x_]:=RightTransversal[x,PermLength[x]];
RightTransversal[GS_GenSet,len_Integer]:=RightTransversal[SchreierSims[{},GS,len],len];
RightTransversal[SGS:StrongGenSet[base_,GS_GenSet],len_Integer]:=Module[{length=len!/OrderOfGroup[SGS],list={},perm,not=NotationOfPerm[SGS]},
While[Length[list]<length,
perm=nosign[First@RightCosetRepresentative[RandomPerm[len,not],len,SGS]];
If[FreeQ[list,perm],AppendTo[list,perm]]
];
PermSort[list]
];
SetNumberOfArguments[RightTransversal,{1,2}];
Protect[RightTransversal];

RightTransversal[TranslatePerm[RiemannSymmetric[{1,2,3,4}],{Perm,4}],4]

RightTransversal[GenSet[Perm[{2,1,3,4}],Perm[{1,2,4,3}]],4]

5.4.2. Left coset transversals

Using the same algorithm we can construct left transversals, where now the canonical representative is chosen according to the order of the inverses.

LeftTransversal[x_]:=LeftTransversal[x,PermLength[x]];
LeftTransversal[GS_GenSet,len_Integer]:=LeftTransversal[SchreierSims[{},GS,len],len];
LeftTransversal[SGS:StrongGenSet[base_,GS_GenSet],len_Integer]:=InversePerm/@RightTransversal[SGS,len];
SetNumberOfArguments[LeftTransversal,{1,2}];
Protect[LeftTransversal];

5.4.3. Double coset transversals

Suppose we are given two subgroups H and K of the group G. Each double coset of H and K in G is a union of several left cosets of K in G and also a union of several right cosets of H in G. There is no overlap among double cosets and hence they also partition the whole group, but they do not all have the same size and therefore there is no Lagrange theorem for double cosets. This makes the computation of the number of double cosets a nontrivial problem.

Here we construct a very simple enumerator (which performs badly when the subgroups H and K have very large index in the group G). As we did with the single coset transversals we shall assume that the group G is the full symmetric group of the degrees of H and K. We also assume that the subgroup K is of the special D-type above because the DoubleCosetRepresentative algorithm is used in the process.

allcases[leftlist_,rightlist_,SGSS_,dummysets_,len_,invert_]:=Module[{list={},lcount,rcount,product,invproduct,canon,frees,slotdummysets},
Do[
product=PermProduct[leftlist[[lcount]],rightlist[[rcount]]];
If[invert,product=InversePerm[product]];
slotdummysets=dummysets/.{Dummyset[m_,pairs_,mQ_]:>DummySet[m,OnPoints[#,product]&/@pairs,mQ],RepeatedSet[list_]:>RepeatedSet[OnPoints[list,product]]};
frees=Complement[Range[len],Union[slotdummysets/.{Dummyset[_,pairs_,_]:>Sequence@@pairs,RepeatedSet[list_]:>list}]];
canon=CanonicalPerm[product,len,SGSS,frees,slotdummysets];
If[invert,canon=InversePerm[canon]];
If[FreeQ[list,canon],AppendTo[list,canon]],
{lcount,1,Length[leftlist]},{rcount,1,Length[rightlist]}];
list]

DoubleTransversal[SGSS_StrongGenSet,dummysets_,len_Integer]:=Module[{transversal,SGSD=SGSOfDummy[dummysets]},
If[OrderOfGroup[SGSS] Length[SGSD[[2]]]<OrderOfGroup[SGSS] Length[SGSS[[2]]],
(* Construct right transversal of S and multiply by D generators on the right *)
transversal=RightTransversal[SGSS,len];
allcases[transversal,SGSD[[2]],SGSS,dummysets,len,True],
(* Construct left transversal of D and multiply by S generators on the left *)
transversal=LeftTransversal[SGSD,len];
allcases[SGSS[[2]],transversal,SGSS,dummysets,len,False]
]
];
DoubleTransversal[Sgroup_StrongGenSet,Dgroup_,n_Integer]:=DoubleTransversal[Sgroup,SGSOfDummySet[Dgroup],n];

6. MathLink connection

(********************** 6. MathLink connection *********************)

{TimeUsed[],SessionTime[],TimeUsed[]/SessionTime[]}

6.0. Comments

If[$ReadingVerbose,Print["Reading section 6: MathLink connection"],Null,Null]

6.1. Check connection

This code installs xperm, printing messages on the process, and sets the global variable $xpermQ to True if the connection was successful, or to False if it failed (perhaps because the system is not Linux , Mac or Windows). Recall that linking is only possible for Windows starting version 6.0 of Mathematica. Actually it is possible to link from previous versions of Mathematica if you have access to the MathLink library ml32i3.dll.

$xActDirectory

xpermConnect:=Module[{version=System`$Version,number=System`$VersionNumber,result,message,dir=StringJoin[$xActDirectory,"/xPerm/mathlink/"]},
Which[

StringMatchQ[version,"*Windows*"],
If[number≥6.0,
message="Connecting to external MinGW executable...";
result=True;
$xpermExecutable=If[StringMatchQ[version,"*64-bit*"],
StringJoin[dir,"xperm.win64"],
StringJoin[dir,"xperm.win32"]
],
message="Pre-6 versions of Mathematica for Windows do not link to MinGW.";
result=False
],

StringMatchQ[version,"*Linux*"],
message="Connecting to external linux executable...";result=True;
$xpermExecutable=If[number≥6.0,
If[StringMatchQ[version,"*64-bit*"],
StringJoin[dir,"xperm.linux.64-bit"],
StringJoin[dir,"xperm.linux.32-bit"]
],
StringJoin[dir,"xperm.linux.pre6"]
],

StringMatchQ[version,"*Mac*"],
message="Connecting to external mac executable...";result=True;$xpermExecutable=StringJoin[$xActDirectory,"/xPerm/mathlink/xperm.mac"],

True,
message="Unsupported system. Contact JMM to compile.";result=False];
Print[message];
If[result,
Check[$xpermLink=Install[$xpermExecutable],result=False];
Print[If[result,"Connection established.","Connection failed."]]
];result
];
$xpermQ=xpermConnect;

We have defined three global variables:

$xpermExecutable

$xpermLink

$xpermQ

LinkPatterns[$xpermLink]

If, for any reason, the connection to the external executable is lost, this function can detect it, reconnect to the executable and repeat the computation. Currently only MLCanonicalPerm is wrapped with it:

SetAttributes[CheckDeadLink,HoldFirst];
CheckDeadLink[expr_]:=Check[expr,Print["Failed computing ",Hold[expr]];$xpermQ=xpermConnect;Print["Reconnected and recomputing"];Module[{result=CheckDeadLink[expr]},Print["Result is ",result];result],LinkObject::linkd];

6.2. Connecting functions

We need to change the formats to connect to xperm.c via MathLink:

6.2.0. FromSign, ToSign

These are two private functions.

We can simulate the existence of signature in permutations of degree n adding two more points n+1 and n+2. Canonical order of those two points means positive signature and reverse order means negative signature. None of the permutations in the problem can mix those two points with the rest of points, so that they form a single orbit.

The full identity will be always a member of the group. We always need to check wheter the reduced identity with negative signature is present or not, using the membership algorithm.

ToSign[perm:Perm[list_],len_]:=If[OnPoints[len+1,perm]===len+1,1,-1]Perm[Drop[list,-2]];
ToSign[perm_Rules,len_]:=If[OnPoints[len+1,perm]===len+1,1,-1]DeleteCases[perm,_[___,len+1,___]];
ToSign[perm_Cycles,len_]:=If[OnPoints[len+1,perm]===len+1,1,-1]DeleteCases[perm,{___,len+1,___}];
ToSign[perm:Images[list_],len_]:=If[OnPoints[len+1,perm]===len+1,1,-1]Images[Drop[list,-2]];
ToSign[GS_GenSet,len_]:=ToSign[#,len]&/@GS;
ToSign[StrongGenSet[base_,GS_],len_]:=StrongGenSet[DeleteCases[base,(len+1)|(len+2)],ToSign[GS,len]];

Note the pattern _[___,len+1,___] instead of the more specific Rule[___,len+1,___]. The latter complains about Rule having three arguments instead of two.

Examples:

ToSign[Perm[{3,4,5,1,2,6,7}],5]

ToSign[Perm[{3,4,5,1,2,7,6}],5]

ToSign[Rules[1->3,3->1],5]

ToSign[Rules[1->3,3->1,6->7,7->6],5]

ToSign[Cycles[{1,3}],5]

ToSign[Cycles[{1,3},{6,7}],5]

ToSign[Images[{1,3,2,4,5,6,7}],5]

ToSign[Images[{1,3,2,4,5,7,6}],5]

Inverse function:

FromSign[s_.Perm[perm_],len_]:=Perm[perm~Join~(len+If[s>0,{1,2},{2,1}])];
FromSign[s_.Rules[rules___],len_]:=Rules@@If[s<0,{rules,len+1->len+2,len+2->len+1},{rules}];
FromSign[s_.Cycles[cycs___],len_]:=Cycles@@If[s<0,{cycs,len+{1,2}},{cycs}];
FromSign[s_.Images[perm_],len_]:=Images[perm~Join~(len+If[s>0,{1,2},{2,1}])];
FromSign[GS_GenSet,len_]:=FromSign[#,len]&/@GS;
FromSign[StrongGenSet[base_,GS_],len_]:=StrongGenSet[If[!FreeQ[GS,-1],Append[base,len+1],base],FromSign[GS,len]];

Examples:

FromSign[-Perm[{3,4,5,1,2}],5]

FromSign[-Rules[1->3,3->1],5]

FromSign[-Cycles[{1,3}],5]

FromSign[-Images[{1,3,2,4,5}],5]

This function converts permutations into lists of images, including the final sign cycle, and generating sets into a single list with all permutations in Images notation with length len+2.

toimagelist[len_][GS_GenSet]:=Flatten@Apply[List,toimagelist[len]/@GS];
toimagelist[len_][sperm_]:=Identity@@TranslatePerm[FromSign[sperm,len],{Images,len+2}];

6.2.1. MathLinkOrbit

Input to MLOrbit:

MathLinkOrbit[point_,GS_GenSet]:=With[{len=Max[point,PermDeg[GS]]},MLOrbit[point,toimagelist[len][GS],len+2]];

6.2.2. MathLinkSchreierSims

Input to MLSchreierSims:

MathLinkSchreierSims[initbase_List,GS_GenSet,len_Integer,options___]:=TranslatePerm[ToSign[MLSchreierSims[initbase,toimagelist[len][GS],len+2],len],NotationOfPerm[GS]];

Output from MLSchreierSims. Note the special role of the Images notation:

Unprotect[StrongGenSet];
StrongGenSet[base:{___Integer},GS:{___Integer},len_Integer]:=(If[Select[base,#>len&]=!={},Print["Computed SGS contains -Cycles[]."]];StrongGenSet[base,GenSet@@(Images/@Partition[GS,len])]);
Protect[StrongGenSet];

StrongGenSet[{1,3},{2,1,3,4,3,4,1,2,1,2,4,3},4]

6.2.3. MathLinkCanonicalPerm

We need to translate the DummySet and RepeatedSet structures into lists of numbers:

TranslateSet[DummySet[_,dums_,sym_]]:={Flatten[{dums}],sym};
TranslateSet[RepeatedSet[list_]]:=Flatten[{list}];
transpose[{}]:={{},{}};
transpose[x_]:=Transpose[x];
TranslateSet[list_List]:=TranslateSetStructure[{transpose[TranslateSet/@Cases[list,DummySet[_,{__},_]]],TranslateSet/@Cases[list,_RepeatedSet]}];
TranslateSetStructure[{{dummysets_,mQ_},repeatedsets_}]:=Sequence[Length/@dummysets,Flatten[dummysets],mQ,Length/@repeatedsets,Flatten[repeatedsets]];

Transpose is a horrible function in cases like this, because it is not idempotent:

Transpose[{{},{}}]

Transpose[{}]

We replace it by transpose:

transpose[{{},{}}]

transpose[%]

Trivial cases:

TranslateSet[{}]

TranslateSet[{DummySet[hello1,{},0]}]

TranslateSet[{DummySet[hello1,{},0],DummySet[hello2,{},0]}]

Nontrivial cases:

TranslateSet[{DummySet[hello1,{},0],DummySet[hello2,{{1,2},{3,4},{5,6}},1]}]

TranslateSet[{DummySet[hello1,{},0],RepeatedSet[{10,11,12,13}],DummySet[hello2,{{1,2},{3,4},{5,6}},1],RepeatedSet[{7,8,9}]}]

TranslateSet[{DummySet[hello1,{},0],RepeatedSet[{10,11,12,13}],DummySet[hello2,{{1,2},{3,4},{5,6}},1],RepeatedSet[{7,8,9}],DummySet[hello3,{{14,15},{16,17}},-1]}]

This function detects whether we are supplying a generating set or a strong generating set, and proceeds accordingly:

tosgslist[GS_GenSet,len_,ob_]:=Sequence[0,If[ob,Range[len],{}],toimagelist[len][GS]];
tosgslist[StrongGenSet[base_List,GS_GenSet],len_,ob_]:=Sequence[1,If[ob,takenext[base,{},Range[len]],Append[base,len+1]],toimagelist[len][GS]];

Input to MLCanonicalPerm. There are two possibilities depending on whether we use a GS or a SGS:

Input structure:
	MathLinkCanonicalPerm[
		permutation: a list in images notation with the two last points encoding the sign
		length: the original length of the signed permutation plus two more points for the sign cycle
		
		SGSswitch: it is 0 when the input is a GS and 1 when the input is a SGS
		base: this is the base of the SGS, or {} when we supply only a GS
		ob: it is 1 when we have reordered the base and 0 otherwise
		GS: this is the given GS or the GS of the SGS, translated into the multiple lists notation
		
		frees: list of free indices (those not participating in the D group)
		
		vds: list of lengths of the nds dummy sets (twice the number of pairs)
		dummies: list of pairs of indices and the symmetry integer (all flattened)
		mQ: list of symmetry integers
		vrs: list of lengths of the nrs dummy sets
		repes: list of repeated indices (all flattened)
	]

Options[CanonicalPerm]

MathLinkCanonicalPerm[sperm_?PermQ,len_Integer,sgs_,frees_List,sets_List,options:OptionsPattern[]]:=Module[{verb,ob,tmp},
{verb,ob}=OptionValue[CanonicalPerm,{options},{xPermVerbose,OrderedBase}];
tmp=TMPHead[
toimagelist[len][sperm],len+2,
tosgslist[sgs,len,ob],
frees,
TranslateSet[sets]
];
If[verb,Print["MathLinkCanonicalPerm:: ",ReplacePart[TMPHead[tmp],{HoldForm,MLCanonicalPerm},{{0},{1,0}},{{1},{2}}]]];
CheckDeadLink[Apply[MLCanonicalPerm,tmp]]/.Images[{(0)..}]->0]

6.2.4. MathLinkSetStabilizer

MathLinkSetStabilizer[list_List,len_Integer,sgs_,options:OptionsPattern[]]:=TranslatePerm[ToSign[
Module[{verb,ob,tmp},
{verb,ob}=OptionValue[SetStabilizer,{options},{xPermVerbose,OrderedBase}];
tmp=Drop[TMPHead[
list,
len+2,
tosgslist[sgs,len,ob]
],{3}];
If[verb,Print["MathLinkSetStabilizer:: ",ReplacePart[TMPHead[tmp],{HoldForm,MathLinkSetStabilizer},{{0},{1,0}},{{1},{2}}]]];
CheckDeadLink[Apply[MLSetStabilizer,tmp]]/.Images[{(0)..}]->0
],
len],NotationOfPerm[sgs]];

6.3. Examples

test[CanonicalPerm,
CanonicalPerm[Images[{2,1,4,3}],4,GenSet[-Images[{3,4,2,1}],Images[{4,3,2,1}]],{1,2,3,4},{DummySet[M2,{},0]},MathLink->False],
Images[{1,2,3,4}]
]//AbsoluteTiming

test[CanonicalPerm,
CanonicalPerm[Images[{2,1,4,3}],4,GenSet[-Images[{3,4,2,1}],Images[{4,3,2,1}]],{1,2,3,4},{DummySet[M2,{},0]},xPermVerbose->True,MathLink->True],
Images[{1,2,3,4}]
]//AbsoluteTiming

test[CanonicalPerm,
CanonicalPerm[Images[{2,1,4,3}],4,GenSet[-Images[{3,4,2,1}],Images[{4,3,2,1}]],{1,2,3,4},{DummySet[M2,{},0]},MathLink->True],
Images[{1,2,3,4}]
]//AbsoluteTiming

test[CanonicalPerm,
CanonicalPerm[Images[{2,1,4,3}],4,StrongGenSet[{1,3},GenSet[-Images[{3,4,2,1}],Images[{4,3,2,1}]]],{1,2,3,4},{DummySet[M2,{},0]},MathLink->False],
Images[{1,2,3,4}]
]//AbsoluteTiming

test[CanonicalPerm,
CanonicalPerm[Images[{2,1,4,3}],4,StrongGenSet[{1,3},GenSet[-Images[{3,4,2,1}],Images[{4,3,2,1}]]],{1,2,3,4},{DummySet[M2,{},0]},MathLink->True],
Images[{1,2,3,4}]
]//AbsoluteTiming

test[CanonicalPerm,
CanonicalPerm[Images[{5,7,3,4,1,2,8,6}],8,StrongGenSet[{1,3,5},GenSet[Cycles[{1,2}],Cycles[{3,4}],Cycles[{5,6}],Cycles[{1,3},{2,4}]]],{1,2},{DummySet[M4,{{5,8},{3,6},{4,7}},1]},MathLink->False],
Images[{3,4,5,6,1,2,7,8}]
]

test[CanonicalPerm,
CanonicalPerm[Images[{5,7,3,4,1,2,8,6}],8,StrongGenSet[{1,3,5},GenSet[Cycles[{1,2}],Cycles[{3,4}],Cycles[{5,6}],Cycles[{1,3},{2,4}]]],{1,2},{DummySet[M4,{{5,8},{3,6},{4,7}},1]},MathLink->True],
Images[{3,4,5,6,1,2,7,8}]
]

Renato's example in paper II:

sperm=TranslatePerm[Perm[{3,2,7,11,4,5,10,6,8,9,12,1}],Images]
GS=TranslatePerm[GenSet[-Perm[{2,1,3,4,5,6,7,8,9,10,11,12}],-Perm[{1,2,4,3,5,6,7,8,9,10,11,12}],-Perm[{1,2,3,4,6,5,7,8,9,10,11,12}],-Perm[{1,2,3,4,5,6,8,7,9,10,11,12}],-Perm[{1,2,3,4,5,6,7,8,10,9,11,12}],-Perm[{1,2,3,4,5,6,7,8,9,10,12,11}],Perm[{3,4,1,2,5,6,7,8,9,10,11,12}],Perm[{1,2,3,4,7,8,5,6,9,10,11,12}],Perm[{1,2,3,4,5,6,7,8,11,12,9,10}],Perm[{1,2,3,4,9,10,11,12,5,6,7,8}]],Images];

Here all indices are free:

test[CanonicalPerm,
CanonicalPerm[sperm,12,GS,Range[1,12],{},TimeVerbose->True,MathLink->False],
Images[{1,5,2,12,3,9,6,8,4,11,7,10}]
]//AbsoluteTiming

test[CanonicalPerm,
CanonicalPerm[sperm,12,GS,Range[1,12],{},MathLink->True],
Images[{1,5,2,12,3,9,6,8,4,11,7,10}]
]//AbsoluteTiming

Here we do not order the base:

test[CanonicalPerm,
CanonicalPerm[sperm,12,GS,Range[1,12],{},TimeVerbose->True,OrderedBase->False,MathLink->False],
-Images[{1,5,2,12,3,9,6,8,7,10,11,4}]
]//AbsoluteTiming

test[CanonicalPerm,
CanonicalPerm[sperm,12,GS,Range[1,12],{},OrderedBase->False,MathLink->True],
-Images[{1,5,2,12,3,9,6,8,7,10,11,4}]
]//AbsoluteTiming

Here we include the dummy indices.

test[CanonicalPerm,
CanonicalPerm[sperm,12,GS,{1,2},{DummySet[M,{{3,4},{5,6},{7,8},{9,10},{11,12}},1]},TimeVerbose->True,MathLink->False],
-Images[{1,3,2,5,4,7,9,11,6,10,8,12}]
]

test[CanonicalPerm,CanonicalPerm[sperm,12,GS,{1,2},{DummySet[M,{{3,4},{5,6},{7,8},{9,10},{11,12}},1]},MathLink->True],
-Images[{1,3,2,5,4,7,9,11,6,10,8,12}]
]//AbsoluteTiming

Here we do not order the base:

test[CanonicalPerm,
CanonicalPerm[sperm,12,GS,{1,2},{DummySet[M,{{3,4},{5,6},{7,8},{9,10},{11,12}},1]},OrderedBase->False,TimeVerbose->True,MathLink->False],
-Images[{1,3,2,6,5,10,7,12,9,8,11,4}]
]

test[CanonicalPerm,
CanonicalPerm[sperm,12,GS,{1,2},{DummySet[M,{{3,4},{5,6},{7,8},{9,10},{11,12}},1]},OrderedBase->False,MathLink->True],
-Images[{1,3,2,6,5,10,7,12,9,8,11,4}]
]//AbsoluteTiming

"Difficult example":

test[CanonicalPerm,
CanonicalPerm[Images[{2,1}],2,GenSet[-Cycles[{1,2}]],{},{DummySet[M2,{{2,1}},1]},MathLink->False],
0
]

test[CanonicalPerm,
CanonicalPerm[Images[{2,1}],2,GenSet[-Cycles[{1,2}]],{},{DummySet[M2,{{2,1}},1]},MathLink->True],
0
]

Now with repeated indices. Suppose we work with a Riemann tensor:

SGS=StrongGenSet[{1,3},GenSet[-Perm[{2,1,3,4}],-Perm[{1,2,4,3}],Perm[{3,4,1,2}]]];

Here the canonical list is {a,b,1,1} for instance. Then the permutation gives R^{1,1,a,b}.

test[CanonicalPerm,
CanonicalPerm[Perm[{3,4,1,2}],4,SGS,{1,2},{RepeatedSet[{3,4}]},MathLink->False],
0
]

test[CanonicalPerm,
CanonicalPerm[Perm[{3,4,1,2}],4,SGS,{1,2},{DummySet[TM,{{3,4}},1]},xPermVerbose->True,MathLink->True],
0
]

test[CanonicalPerm,
CanonicalPerm[Perm[{3,4,1,2}],4,SGS,{1,2},{RepeatedSet[{3,4}]},xPermVerbose->True,OrderedBase->False,MathLink->True],
0
]

TranslateSet[{RepeatedSet[{3,4}]}]

Simple example with SchreierSims, including negative permutations:

SchreierSims[{},GenSet[-Cycles[{1,2}],Cycles[{1,3},{2,4}]],4,MathLink->False]

OrderOfGroup[%,MathLink->True]

SchreierSims[{},GenSet[-Cycles[{1,2}],Cycles[{1,3},{2,4}]],4,MathLink->True]

OrderOfGroup[%,MathLink->True]

7. End

7.1. End private and package

End[]

EndPackage[]

On[System`Cycles::shdw];

?xAct`xPerm`*

?xAct`xPerm`Private`*

Time spent in executing the notebook:

SessionTime[]

TimeUsed[]

Checks:

xAct`xPerm`Private`alltests[]

{Length@%,And@@%}

Flatten@Position[%%,False]

DownValues[xAct`xPerm`Private`testresult][[%]]