Using GAP Effectively: Some Tips and Pitfalls

Based on Steve Linton's talk at the Second CoDiMa Training School in Computational Discrete Mathematics (ICMS, Edinburgh, 17-21 October 2016, https://www.codima.ac.uk/school2016/)

Philosophy

  • GAP has two somewhat contradictory design goals
    • to allow users to pose questions in a way that seems natural to a working mathematician and get answers
    • to allow the expert computational mathematician to implement and apply the most advanced techniques to solve hard problems
  • The first is achieved to a limited extent.

Motivating problem: find an element of $S_9$ which is NOT an involution

In [1]:
Filtered(Elements(SymmetricGroup(9)), x-> x*x <> ())[1];
Out[1]:
(7,8,9)
  • Replace 9 by 10 and it is already visibly longer
In [2]:
Filtered(Elements(SymmetricGroup(10)), x-> x*x <> ())[1];
Out[2]:
(8,9,10)
  • Put n:=15 and you quickly run out of memory: $15!$ is roughly $1.3 \times 10^{12}$
  • Of course GAP is capable of doing the latter calculation
  • But you need to acquire certain expertise to ask questions in the right way
  • This talk is about how to start thinking like an expert

If you don't need it, then don't store it!

n:=9;; Filtered( Elements(SymmetricGroup(n) ),
                 x-> x*x <> () )[1];
  • First this computes and stores the full list of elements of $S_n$
  • Then it checks each of them to see if it has order 2
  • It stores a second list of all of those which don’t
  • Finally it returns the first one
  • Instead, we can stop looking when we find one. GAP even provides a built in function to do this:
In [3]:
First(Elements(SymmetricGroup(9)), x-> x*x <> ());
Out[3]:
(7,8,9)
  • Stopping things as soon as possible is an important principle
  • In this case, however, the real problem is computing and storing all the elements of $S_n$
  • Let's explore some other alternatives

Enumerators

  • Enumerator returns a list of elements of a domain
    • this list may be virtual
    • also EnumeratorSorted — but only if you need it
  • For many objects it is quick to construct, but may be slower to access
In [4]:
# Do not remove the 2nd semicolon below yet
# See https://github.com/gap-packages/JupyterKernel/issues/72
e := Enumerator(SymmetricGroup(99));;
In [5]:
Length(e);
Out[5]:
933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000
In [6]:
e[10^100];
Out[6]:
(2,60,99,55,54,65,7,16,18,32,70,15,5,37,43,97,19,31,66,30,90,17,29,85,28,67,27,62,26,34,52,59)(3,44,73,47,95,45,51,68,50,86,49,83,40,36,81,35,93,12,76,11,75,10,46,9,96,8,53,42,41,22,78,21,38,20,24,63,23,48,39,56,4,6,58,14,80,13,25,33)

See also EnumeratorOfTuples and EnumeratorOfCombinations (both documented) and EnumeratorOfCartesianProduct (currently undocumented)

Iterators

  • Even an Enumerator can be too heavyweight
    • sometimes you don’t need to even number the elements, or know how many there are
  • For this GAP has Iterators
    • IsDoneIterator and NextIterator operations
In [8]:
n := 9;; i := Iterator(SymmetricGroup(n));;
In [9]:
while not IsDoneIterator(i) do 
    x := NextIterator(i); 
    if x*x <> () then 
        break; 
    fi; 
od;
In [10]:
x;
Out[10]:
(7,9,8)

Or more concisely, thanks to some built-in magic:

In [12]:
n := 9;; 
for x in SymmetricGroup(n) do 
    if x*x <> () then 
        break; 
    fi; 
od;
In [13]:
x;
Out[13]:
(7,9,8)

Or even shorter:

In [15]:
n := 9;; First(SymmetricGroup(n), x->x*x <> ());
Out[15]:
(7,9,8)

Randomness

  • Sometimes you can’t even make an iterator for your group easily, but you know the elements you want exist and are not too rare
  • So make pseudo-random elements of the group until you find one
In [16]:
g := SL(10,3);
Out[16]:
SL(10,3)
In [17]:
repeat 
    x := PseudoRandom(g); 
until Order(x) = (3^10-1)/2;
In [18]:
x;
Out[18]:
[ [ 0*Z(3), Z(3)^0, Z(3), Z(3), Z(3)^0, Z(3)^0, 0*Z(3), Z(3)^0, Z(3), Z(3)^0 ], [ Z(3), Z(3), Z(3), Z(3)^0, 0*Z(3), Z(3)^0, Z(3)^0, Z(3), Z(3)^0, 0*Z(3) ], [ 0*Z(3), Z(3)^0, 0*Z(3), 0*Z(3), Z(3), 0*Z(3), Z(3), 0*Z(3), Z(3), Z(3) ], [ Z(3), 0*Z(3), Z(3), Z(3)^0, Z(3)^0, Z(3), Z(3)^0, Z(3), Z(3)^0, 0*Z(3) ], [ Z(3)^0, Z(3)^0, Z(3)^0, Z(3), 0*Z(3), 0*Z(3), 0*Z(3), 0*Z(3), 0*Z(3), Z(3) ], [ Z(3), Z(3)^0, 0*Z(3), Z(3)^0, Z(3), 0*Z(3), Z(3)^0, Z(3), 0*Z(3), Z(3)^0 ], [ Z(3), Z(3)^0, Z(3)^0, Z(3)^0, Z(3), Z(3)^0, 0*Z(3), 0*Z(3), Z(3), Z(3)^0 ], [ Z(3)^0, Z(3)^0, Z(3), Z(3), Z(3)^0, Z(3), 0*Z(3), Z(3), Z(3), Z(3) ], [ Z(3)^0, Z(3)^0, 0*Z(3), Z(3), Z(3)^0, Z(3), Z(3), Z(3), Z(3), Z(3)^0 ], [ 0*Z(3), Z(3), 0*Z(3), Z(3), Z(3), Z(3)^0, Z(3), Z(3), Z(3)^0, Z(3) ] ]
In [19]:
Display(x);
 . 1 2 2 1 1 . 1 2 1
 2 2 2 1 . 1 1 2 1 .
 . 1 . . 2 . 2 . 2 2
 2 . 2 1 1 2 1 2 1 .
 1 1 1 2 . . . . . 2
 2 1 . 1 2 . 1 2 . 1
 2 1 1 1 2 1 . . 2 1
 1 1 2 2 1 2 . 2 2 2
 1 1 . 2 1 2 2 2 2 1
 . 2 . 2 2 1 2 2 1 2

But is searching through all the elements the right thing to do in the first place?

  • Element order is a conjugacy invariant
  • For many groups there are ways of finding conjugacy class representatives that are faster than listing all elements
    • or they might be already known and stored
In [21]:
n := 9;; 
Representative(First(ConjugacyClasses(SymmetricGroup(n)),
    c->Representative(c)^2 <> ()));
Out[21]:
(1,2,3)
  • This is one of the most powerful techniques, especially for non-abelian simple groups and things close to them
  • Of course if you are really working in $S_n$ you can simply construct the answer as a permutation
In [22]:
First(SymmetricGroup(11), 
    x-> OnTuples([1,2,3,4,5],x) = [1,3,5,7,9] and 
        Order(x) = 7);
Out[22]:
(2,3,5,9,4,7,11)
  • For values larger than 12, this gets slow
    • because it searches lots of elements that fix 2 before it looks at anything that moves 1 to 2
  • Use a bit of maths
    • the elements that maps [1,2,3,4,5] to [1,3,5,7,9] lie in a coset of a sequence stabilizer
In [24]:
g := SymmetricGroup(12);; s := Stabilizer(g,[1,2,3,4,5],OnTuples);;
In [25]:
r := RepresentativeAction(g,[1,2,3,4,5],[1,3,5,7,9],OnTuples);
Out[25]:
(2,3,5,9,4,7)
In [26]:
First(s,x->Order(x*r) = 7)*r;
Out[26]:
(2,3,5,9,12,4,7)

General Principles

Searching for an element in a group

  • Don’t write down the list of elements first
  • Stop when you’ve found it
  • Stop looking at other elements as soon as you know they’re not it
    • order of a matrix can be large and a bit slow to compute
    • if all you care about is whether it is 2, just check for IsOne(x*x) and not IsOne(x)

Searching for an element in a group (2)

  • Try to identify a subgroup, or coset or conjugacy class that it lies in
    • remember Sylow subgroups!
    • automorphism group sometimes helps too
  • Search only in there

Seaching for a subgroup

  • Even worse — quite small groups can have very many subgroups
  • Some kinds that are eas(ier) to find
    • Cyclic subgroups (via ConjugacyClasses)
    • NormalSubgroups
    • Derived, Lower Central etc. series
    • Sylow subgroups
    • ...

Seaching for a subgroup (2)

  • Some kinds that are eas(ier) to find
    • ...
    • Maximal subgroups (for some groups)
    • MaximalSubgroups will return all subgroups. You are likely to want ony MaximalSubgroupClassReps
  • Ask yourself if one of these lists might include the one you want, or at least help you on your way

Searching for multiple elements

Conjecture: $U_3(3)$ cannot be generated by three involutions

In [27]:
g := PSU(3,3);
Out[27]:
<permutation group of size 6048 with 2 generators>

So we know some things not to do

  • list all 216G triples of elements of $U_3(3)$ and filter out all the ones that generate the group and consist of involutions
  • use IteratorOfTuples to run through all 216G ...
  • use IteratorOfCombinations to run through 36G of unordered triples
  • the same, but test for involutions first
    • would take a few hours on a laptop

What to do

  • find the involutions first (there are just 63 of them) and run over triples
In [28]:
invs := Filtered(g, x -> IsOne(x*x) and not IsOne(x));;
In [29]:
Length(invs);
Out[29]:
63
In [32]:
iter := IteratorOfCombinations(invs,3);; ct := 0;;
while not IsDoneIterator(iter) do
    x := NextIterator(iter);
    if Subgroup(g,x) = g then 
        break; 
    fi; 
    ct := ct+1; 
od;
In [33]:
ct;
Out[33]:
39711
In [34]:
Binomial(63,3);
Out[34]:
39711

Searching for multiple elements

  • We still haven’t used conjugacy
  • Could choose the 1st involution to be a conjugacy class representative
    • there is only one conjugacy class of involutions
    • reduce search from Binomial(63,3) to Binomial(62,2)
  • The 2nd involution can be chosen up to conjugacy in the centraliser of the first one
    • just 4 cases to consider, so the search is now $4 \times 61$ cases
  • Of course the 3rd one can be chosen up to conjugacy in the normaliser of the subgroup generated by the first two
  • If the things you are searching for are not all the same, then the order in which you look at them also matters

Morpheus

  • This type of search for sequences of elements that generate something is nicely implemented by Alexander Hulpke in a part of the GAP library called Morpheus
  • There are various functions that access morpheus documented in the library under "Searching for Homomorphisms"
  • Our example is asking whether $U_3(3)$ is a quotient of the free product of three cyclic groups of order 2
In [35]:
g:=PSU(3,3);;
In [36]:
F:=FreeGroup(3);
Out[36]:
<group with 3 generators>
In [37]:
F:=F/[F.1^2,F.2^2,F.3^2];
Out[37]:
<group with 3 generators>
In [38]:
GQuotients(F,g);
Out[38]:
[  ]

Another Morpheus example

In [39]:
F:=FreeGroup(2);
Out[39]:
<group with 2 generators>
In [40]:
F:=F/[F.1^2,F.2^6];
Out[40]:
<group with 2 generators>
In [41]:
qs:=GQuotients(F,g);
Out[41]:
[ [ f1, f2 ] -> [ (2,4)(5,9)(6,7)(8,10)(11,21)(12,22)(13,20)(14,24)(15,28)(16,27)(17,25)(18,26)(19,23)(38,60)(39,64)(40,59)(41,61)(42,63)(43,58)(44,57)(45,56)(46,62)(47,73)(48,68)(49,69)(50,72)(51,71)(52,65)(53,67)(54,66)(55,70)(74,89)(75,88)(76,90)(77,83)(78,84)(79,87)(80,86)(81,91)(82,85), (1,43,21,15,26,58)(2,91,73)(3,22,89,24,66,25)(4,33)(5,47,13,87,36,12)(6,57,45,78,76,83)(7,17,29,69,16,81)(8,10,77,60,42,50)(9,72,52,51,56,44)(11,48,80,18,67,88)(19,64,46)(20,74,23,54,27,32)(30,59,40,85,86,75)(31,90,62,39,65,34)(35,53,68,70,61,41)(37,82,55)(38,63)(49,84)(71,79) ], [ f1, f2 ] -> [ (2,5)(3,10)(4,6)(7,9)(11,85)(12,83)(13,84)(14,91)(15,86)(16,89)(17,90)(18,88)(19,87)(20,54)(21,53)(22,52)(23,49)(24,47)(25,50)(26,55)(27,51)(28,48)(29,46)(30,41)(31,42)(32,45)(33,44)(34,38)(35,40)(36,39)(37,43)(65,77)(66,78)(67,82)(68,80)(69,79)(70,75)(71,74)(72,76)(73,81), (1,43,21,15,26,58)(2,91,73)(3,22,89,24,66,25)(4,33)(5,47,13,87,36,12)(6,57,45,78,76,83)(7,17,29,69,16,81)(8,10,77,60,42,50)(9,72,52,51,56,44)(11,48,80,18,67,88)(19,64,46)(20,74,23,54,27,32)(30,59,40,85,86,75)(31,90,62,39,65,34)(35,53,68,70,61,41)(37,82,55)(38,63)(49,84)(71,79) ] ]
In [42]:
Length(qs);
Out[42]:
2
  • So $U_3(3)$ is (2,6)-generated in two distinct ways
  • Presented as homomorphisms — easy to recover the generators if you want them
  • Other Morpheus functions: AllHomomorphisms, AutomorphismGroup, IsomorphicSubgroups
  • A powerful tool for many purposes

Working in the right Group

  • Mathematicians are very sloppy: they constantly identify isomorphic groups
  • So $A_5$ "is" $PSL(2,5)$ and $SL(2,4)$ and $\langle a,b \mid a^2=b^3=(ab)^5=1 \rangle$ and $\langle (1,3,6,2,4), (1,2,3)(4,5,6) \rangle$
  • Computationally these are different, so you must choose the right one to work in
  • Two tools for moving between them: homomorphisms and straight-line programs

Finitely Presented Groups

  • Lots of functionality in GAP for fp groups — mostly to do with identifying unknown ones
  • Lots of textbooks that define groups by presentations, e.g. $$D_{2n} = \langle a,b \mid a^n = (ab)^2 = b^2 = 1 \rangle $$
  • GAP supports some general group theoretic computation with fp groups that turn out to be finite
  • But it’s usually the wrong way to do things
In [43]:
f := FreeGroup("a","b");
Out[43]:
<group with 2 generators>
In [44]:
AssignGeneratorVariables(f);
#I  Assigned the global variables [ a, b ]
In [45]:
g := f/[a^2,b^3,(a*b)^7, Comm(a,b)^8];
Out[45]:
<group with 2 generators>
In [46]:
Sum(Elements(g), Order);
Out[46]:
66655
In [47]:
x := Random(g);
Out[47]:
<object>
In [48]:
Print(x);
(b*a^-1*b^-1*a^-1)^3*b*a^-1*b^-1*(a^-1*b^-1*(a^-1*b^-1*a^-1*b)^3*a^-1*b^-1*a^-\
1)^2*b^-1*a^-1*(a^-1*b^-1)^2*a^-1
In [49]:
g := f/[a^2,b^3,(a*b)^7, Comm(a,b)^8];
Out[49]:
<group with 2 generators>
In [50]:
phi := IsomorphismPermGroup(g);
Out[50]:
<object>
In [51]:
Print(phi);
GroupHomomorphismByImages( Group( [ a, b ] ), Group( 
[ ( 1, 2)( 3, 5)( 4, 6)( 7,11)( 8,12)( 9,13)(10,14)(16,20)(17,21)(18,22)
    (19,23)(25,29)(26,27)(28,30)(31,34)(32,35)(33,36)(37,41)(38,42)(39,43)
    (40,44)(45,50)(46,51)(48,52)(49,53)(54,56), 
  ( 2, 3, 4)( 5, 7, 8)( 6, 9,10)(11,15,14)(12,16,17)(13,18,19)(20,24,23)
    (21,25,26)(22,27,28)(29,31,32)(30,33,34)(35,37,38)(36,39,40)(41,45,46)
    (42,47,43)(44,48,49)(50,53,54)(51,55,52) ] ), [ a, b ], 
[ ( 1, 2)( 3, 5)( 4, 6)( 7,11)( 8,12)( 9,13)(10,14)(16,20)(17,21)(18,22)
    (19,23)(25,29)(26,27)(28,30)(31,34)(32,35)(33,36)(37,41)(38,42)(39,43)
    (40,44)(45,50)(46,51)(48,52)(49,53)(54,56), 
  ( 2, 3, 4)( 5, 7, 8)( 6, 9,10)(11,15,14)(12,16,17)(13,18,19)(20,24,23)
    (21,25,26)(22,27,28)(29,31,32)(30,33,34)(35,37,38)(36,39,40)(41,45,46)
    (42,47,43)(44,48,49)(50,53,54)(51,55,52) ] )
In [52]:
h := ImagesSource(phi);
Out[52]:
<permutation group of size 10752 with 2 generators>
In [53]:
Sum(Elements(h), Order);
Out[53]:
66655
In [54]:
x := Random(h);
Out[54]:
(1,9,31)(2,5,39)(3,27,20)(4,35,55)(6,21,25)(7,37,30)(8,24,51)(11,41,47)(12,28,46)(13,43,38)(14,18,22)(15,53,17)(16,40,23)(19,54,45)(26,56,44)(29,49,32)(33,52,36)(34,42,50)
In [55]:
y := PreImagesRepresentative(phi,x);
Out[55]:
<object>
In [56]:
Print(y);
b^-1*a^-1*b*((a^-1*b^-1*a^-1*b)^3*a^-1*b^-1)^2*a^-1*(b^-1*(a*b)^2*a)^2*(b^-1*a\
)^2*b

Other Isomorphism Constructors

  • Isomorphism[Special]PcGroup
    • pcgroups are usually the fastest representation for solvable groups
  • IsomorphismFpGroup
    • basically only if you want a presentation of your group
  • SmallerDegreePermRep
    • heuristic
  • GAP will sometimes do this for you
    • see ?NiceMonomorphism or ?NiceObject
    • but it can be better to do it by hand

A Few Homomorphism Operations

  • Part of general mapping (relation) machinery
  • Source and Range (domain and codomain)
    • given when the morphism is constructed
    • morphism does not need to be total or onto, so they may be bigger than you expect
    • ImagesSource and PreImagesRange may be what you want
  • Image specialised to ImageElm and ImagesSet
    • which don’t check that the input is in the source
  • PreImagesRepresentative gives just ONE preimage
  • InverseGeneralMapping
  • CompositionMapping
In [57]:
g:=Group((1,2,3,4),(1,2),(5,6,7));
Out[57]:
Group([ (1,2,3,4), (1,2), (5,6,7) ])
In [58]:
iso:=IsomorphismPcGroup(g);
Out[58]:
<object>
In [59]:
h:=Image(iso);
Out[59]:
<group of size 72 with 5 generators>
In [60]:
z:=Centre(h);
Out[60]:
<group with 1 generators>
In [61]:
SetCentre(g,PreImage(iso,z));
In [62]:
cl:=ConjugacyClasses(h);
Out[62]:
[ <object>, <object>, <object>, <object>, <object>, <object>, <object>, <object>, <object>, <object>, <object>, <object>, <object>, <object>, <object> ]
In [63]:
ncl:=[];
Out[63]:
[  ]
In [64]:
for c in cl do
    nc:=ConjugacyClass(g,PreImage(iso,Representative(c)));
    SetSize(nc,Size(c));
    SetStabilizerOfExternalSet(nc, PreImage(iso,StabilizerOfExternalSet(c)));
    Add(ncl,nc);
od;
In [65]:
List(ncl,Size);
Out[65]:
[ 1, 1, 6, 8, 3, 1, 6, 8, 3, 6, 6, 8, 3, 6, 6 ]
In [66]:
SetConjugacyClasses(g,ncl);

Homorphisms in General

  • Even if you can’t find an isomorphism to a nicer group, you may be able to find a homomorphism
  • Solve your problem in the image first and refine
In [67]:
g := Group((1,2),(3,4),(5,6),(7,8),(9,10,11),(11,12,13));
Out[67]:
Group([ (1,2), (3,4), (5,6), (7,8), (9,10,11), (11,12,13) ])
In [69]:
Number(g, x-> Order(x) mod 2 = 1); Size(g);
Out[69]:
45
Out[69]:
960
In [70]:
Orbits(g,MovedPoints(g));
Out[70]:
[ [ 1, 2 ], [ 3, 4 ], [ 5, 6 ], [ 7, 8 ], [ 9, 10, 11, 12, 13 ] ]
In [71]:
phi := ActionHomomorphism(g,[1..8]);
Out[71]:
<object>
In [72]:
Print(phi);
<action homomorphism>
In [73]:
h := ImagesSource(phi);
Out[73]:
Group([ (1,2), (3,4), (5,6), (7,8) ])
In [74]:
odds := Filtered(h, x->Order(x) mod 2 = 1);
Out[74]:
[ () ]
In [75]:
p := PreImagesSet(phi,odds);
Out[75]:
[ (), (11,12,13), (11,13,12), (10,11)(12,13), (10,11,12), (10,11,13), (10,12,11), (10,12,13), (10,12)(11,13), (10,13,11), (10,13,12), (10,13)(11,12), (9,10)(12,13), (9,10)(11,12), (9,10)(11,13), (9,10,11), (9,10,11,12,13), (9,10,11,13,12), (9,10,12,13,11), (9,10,12), (9,10,12,11,13), (9,10,13,12,11), (9,10,13), (9,10,13,11,12), (9,11,10), (9,11,12,13,10), (9,11,13,12,10), (9,11)(12,13), (9,11,12), (9,11,13), (9,11)(10,12), (9,11,10,12,13), (9,11,13,10,12), (9,11)(10,13), (9,11,10,13,12), (9,11,12,10,13), (9,12,13,11,10), (9,12,10), (9,12,11,13,10), (9,12,11), (9,12,13), (9,12)(11,13), (9,12,13,10,11), (9,12)(10,11), (9,12,10,11,13), (9,12,10,13,11), (9,12,11,10,13), (9,12)(10,13), (9,13,12,11,10), (9,13,10), (9,13,11,12,10), (9,13,11), (9,13,12), (9,13)(11,12), (9,13,12,10,11), (9,13)(10,11), (9,13,10,11,12), (9,13,10,12,11), (9,13,11,10,12), (9,13)(10,12) ]
In [77]:
Length(p); Number(p, x->Order(x) mod 2 = 1);
Out[77]:
60
Out[77]:
45

Not all homomorphisms are equal

  • If you just make a GroupHomorphismByImages (by giving images of generators)
    • it can be slow to make because it checks (use GroupHomorphismByImagesNC if you are sure you are right)
    • Image and preimage computation can be slow, or preimages can be “nasty” (long words in FP group)
      • essential because factorisation in terms of generators is not always easy
  • ActionHomorphism is usually good
  • So are most things produced by IsomorphismXXXGroup

Random Tips 1

Avoid long lists of mutable objects

  • since the objects in the list might change “under its feet” the list can’t remember
    • whether it’s sorted
    • whether the entries are all from the same family
  • so whenever you try and search it or call an operation on it, it has to look at every element
    • can become very slow
  • lists of immutable objects are much better
  • sorted lists of immutable comparable objects can use binary search

Random Tips 2

  • There are space and time efficient representations of vectors and matrices over finite fields
    • up to order 256 in the kernel
    • bigger fields in package cvec
  • Vectors and matrices are not always in these representations by default
    • among other reasons because deciding whether this vector is “really” over GF(3) or GF(9) requires prescience
  • ConvertToVectorRep(v, q) and ConvertToMatrixRep(m,q) convert in place
  • cvec has its own functions
  • working with large uncompressed vectors or matrices is a bad idea

Further Reading

  • A lot of this talk was taken from Alexander Hulpke’s talk “Using GAP”, especially section 4
  • You can read the original paper by Alexander Hulpke at http://www.math.colostate.edu/~hulpke/paper/gap4tut.pdf
  • A lot of similar ideas are found in Steve Linton's paper "The Art and Science of Computing in Large Groups" (in Bosma & van der Poorten: Computational Algebra and Number Theory, 1995, Springer)