How do I construct … in GAP?

This page is devoted to answering some basic questions along the line
“How do I construct … in GAP?” You may view the html source code
for the GAP commands without the output or GAP prompt.

Please send suggestions, additions, corrections to
David Joyner.

This page itself is under construction…


Questions

How
do I construct a … group?

permutation
dihedral 
cyclicconjugacy classes of a
finitely presented

How
do I … a polynomial?
How do I find the … of a group
representation?
How
do I compute an mod m, where A is …?
Given
a group G, how do I compute … ?

Answers


    • permutation:
      To construct a permutation group, write down generators in disjoint cycle notation,
      put them in a list (i.e., surround them by square brackets), andThe permutation group G generated by the cycles
      (1,2)(3,4) and (1,2,3):
gap> G:=Group((1,2)(3,4),(1,2,3));

Group([ (1,2)(3,4), (1,2,3) ])

This is of course a subgroup of the symmetric group S4 on 4
letters.
Indeed, this G is in fact the alternating group
on four letters, A4.

By virtue of the fact that the permutations generating G employ
integers less than or equal to 4, this group G
is a subgroup of the symmetric group S4 on 4
letters. Some permutation groups have special constructions:

gap> S4:=SymmetricGroup(4);
Sym( [ 1 .. 4 ] )
gap> A4:=AlternatingGroup(4);
Alt( [ 1 .. 4 ] )
gap> IsSubgroup(S4,G);
true
gap> IsSubgroup(A4,G);
true
gap> S3:=SymmetricGroup(3);
Sym( [ 1 .. 3 ] )
gap> IsSubgroup(S3,G);
false


    • dihedral
      To construct a dihedral group, use the special “DihedralGroup” command:
gap> G:=DihedralGroup(6);

gap> Size(G);
6
gap> f:=GeneratorsOfGroup( G );
[ f1, f2 ]
gap> f[1]^2; f[2]^3;
identity of ...
identity of ...
gap> f[1]^2= f[2]^3;
true


  • cyclic group
    To construct a cyclic group, you may
    construct integers mod n:

    gap> R:=ZmodnZ( 12);
    (Integers mod 12)
    gap> a:=Random(R);
    ZmodnZObj( 11, 12 )
    gap> 4*a;
    ZmodnZObj( 8, 12 )
    gap> b:=Random(R);
    ZmodnZObj( 9, 12 )
    gap> a+b;
    ZmodnZObj( 8, 12 )
    

    or use the special “CyclicGroup” command

    gap> G:=CyclicGroup(12);
    pc group of size 12 with 3 generators
    gap> a:=Random(G);
    f3^2
    gap> f:=GeneratorsOfGroup( G );
    [ f1, f2, f3 ]
    gap> f[1]^4;
    f3
    gap> f[1]^12;
    identity of ...
    
    


  • conjugacy:
    The conjugacy classes of a group G are computed using
    the “ConjugacyClasses” command. This is a list
    of classes{x^-1*g*x | x in G}.

    gap> G:=SL(2,7);
    SL(2,7)
    gap> CG:=ConjugacyClasses(G);
    [ [ [ Z(7)^0, 0*Z(7) ], [ 0*Z(7), Z(7)^0 ] ]^G,
      [ [ 0*Z(7), Z(7)^3 ], [ Z(7)^0, Z(7)^5 ] ]^G,
      [ [ 0*Z(7), Z(7)^4 ], [ Z(7)^5, Z(7)^5 ] ]^G,
      [ [ Z(7)^3, 0*Z(7) ], [ 0*Z(7), Z(7)^3 ] ]^G,
      [ [ 0*Z(7), Z(7)^3 ], [ Z(7)^0, Z(7)^2 ] ]^G,
      [ [ 0*Z(7), Z(7)^4 ], [ Z(7)^5, Z(7)^2 ] ]^G,
      [ [ 0*Z(7), Z(7)^3 ], [ Z(7)^0, 0*Z(7) ] ]^G,
      [ [ 0*Z(7), Z(7)^3 ], [ Z(7)^0, Z(7)^4 ] ]^G,
      [ [ 0*Z(7), Z(7)^3 ], [ Z(7)^0, Z(7) ] ]^G,
      [ [ Z(7)^4, 0*Z(7) ], [ 0*Z(7), Z(7)^2 ] ]^G,
      [ [ Z(7)^5, 0*Z(7) ], [ 0*Z(7), Z(7) ] ]^G ]
    gap> g:=Representative(CG[3]); Order(g);
    [ [ 0*Z(7), Z(7)^4 ], [ Z(7)^5, Z(7)^5 ] ]
    14
    gap> g:=Representative(CG[4]); Order(g);
    [ [ Z(7)^3, 0*Z(7) ], [ 0*Z(7), Z(7)^3 ] ]
    2
    gap> g:=Representative(CG[5]); Order(g);
    [ [ 0*Z(7), Z(7)^3 ], [ Z(7)^0, Z(7)^2 ] ]
    7
    gap> g:=Representative(CG[6]); Order(g);
    [ [ 0*Z(7), Z(7)^4 ], [ Z(7)^5, Z(7)^2 ] ]
    7
    gap>                            
    


  • presented
    To construct a finitely presented group in GAP, use the
    “FreeGroup” and “” commands. Here is one example.

    gap> M12 := MathieuGroup( 12 );
    Group([ (1,2,3,4,5,6,7,8,9,10,11), (3,7,11,8)(4,10,5,6), (1,12)(2,11)(3,6)(4,8)(5,9)(7,10) ])
    gap> F := FreeGroup( "a", "b", "c" );
    free group on the generators [ a, b, c ]
    gap> words := [ F.1, F.2 ];
    [ a, b ]
    gap> P := PresentationViaCosetTable( M12, F, words );
    presentation with 3 gens and 10 rels of total length 97
    gap> TzPrintRelators( P );
    #I  1. c^2
    #I  2. b^4
    #I  3. a*c*a*c*a*c
    #I  4. a*b^2*a*b^-2*a*b^-2
    #I  5. a^11
    #I  6. a^2*b*a^-2*b^2*a*b^-1*a^2*b^-1
    #I  7. a*b*a^-1*b*a^-1*b^-1*a*b*a^-1*b*a^-1*b^-1
    #I  8. a^2*b*a^2*b^2*a^-1*b*a^-1*b^-1*a^-1*b^-1
    #I  9. a*b*a*b*a^2*b^-1*a^-1*b^-1*a*c*b*c
    #I  10. a^4*b*a^2*b*a^-2*c*a*b*a^-1*c
    gap> G := FpGroupPresentation( P );
    fp group on the generators [ a, b, c ]
    gap> RelatorsOfFpGroup( G );  
    [ c^2, b^4, a*c*a*c*a*c, a*b^-2*a*b^-2*a*b^-2, a^11, a^2*b*a^-2*b^-2*a*b^-1*a^2*b^-1, a*b*a^-1*b*a^-1*b^-1*a*b*a^-1*b*a^-1*b^-1,
      a^2*b*a^2*b^-2*a^-1*b*a^-1*b^-1*a^-1*b^-1, a*b*a*b*a^2*b^-1*a^-1*b^-1*a*c*b*c, a^4*b*a^2*b*a^-2*c*a*b*a^-1*c ]
    gap> Size(M12);
    95040
    gap> Size(G);
    95040
    gap> IsomorphismGroups(G,M12); 
    ????????
    

    The last command is computationally intensive and requires more
    than the default memory allocation of 256M of RAM.

    Here is another example.

    gap> F := FreeGroup( "a", "b");
    free group on the generators [ a, b ]
    gap> G:=F/[F.1^2,F.2^3,F.1*F.2*F.1^(-1)*F.2^(-1)];
    fp group on the generators [ a, b ]
    gap> Size(G);
    6
    
    


  • rref
    The key command for row reduction is “TriangulizeMat”.
    The following example illustrates the syntax.

    gap> M:=[[1,2,3,4,5],[1,2,1,2,1],[1,1,0,0,0]];
    [ [ 1, 2, 3, 4, 5 ], [ 1, 2, 1, 2, 1 ], [ 1, 1, 0, 0, 0 ] ]
    gap> TriangulizeMat(M);
    gap> M;
    [ [ 1, 0, 0, -1, 1 ], [ 0, 1, 0, 1, -1 ], [ 0, 0, 1, 1, 2 ] ]
    gap> Display(M);
    [ [   1,   0,   0,  -1,   1 ],
      [   0,   1,   0,   1,  -1 ],
      [   0,   0,   1,   1,   2 ] ]
    gap> M:=Z(3)^0*[[1,2,3,4,5],[1,2,1,2,1],[1,1,0,0,0]];
    [ [ Z(3)^0, Z(3), 0*Z(3), Z(3)^0, Z(3) ],
      [ Z(3)^0, Z(3), Z(3)^0, Z(3), Z(3)^0 ],
      [ Z(3)^0, Z(3)^0, 0*Z(3), 0*Z(3), 0*Z(3) ] ]
    gap> TriangulizeMat(M);
    gap> Display(M);
     1 . . 2 1
     . 1 . 1 2
     . . 1 1 2
    gap>
    


  • kernel:
    There are different methods for matrices over the integers and
    matrices over a field.For integer entries, related commands include
    “NullspaceIntMat” and “SolutionNullspaceIntMat”
    in section

    25.1 “Linear equations over the integers and Integral Matrices”

    of the reference manual.

    gap> M:=[[1,2,3],[4,5,6],[7,8,9]];
    [ [ 1, 2, 3 ], [ 4, 5, 6 ], [ 7, 8, 9 ] ]
    gap> NullspaceIntMat(M);
    [ [ 1, -2, 1 ] ]
    gap> SolutionNullspaceIntMat(M,[0,0,1]);
    [ fail, [ [ 1, -2, 1 ] ] ]
    gap> SolutionNullspaceIntMat(M,[0,0,0]);
    [ [ 0, 0, 0 ], [ [ 1, -2, 1 ] ] ]
    gap> SolutionNullspaceIntMat(M,[1,2,3]);
    [ [ 1, 0, 0 ], [ [ 1, -2, 1 ] ] ]
    
    

    Here (0,0,1) is not in the image of M
    (under v-> v*M) but (0,0,0) and (1,2,3) are.

    For field entries, related commands include
    “NullspaceMat” and “TriangulizedNullspaceMat”
    in section

    24.6 “Matrices Representing Linear Equations and the Gaussian Algorithm”

    of the reference manual.

    gap> M:=[[1,2,3],[4,5,6],[7,8,9]];
    [ [ 1, 2, 3 ], [ 4, 5, 6 ], [ 7, 8, 9 ] ]
    gap> NullspaceMat(M);
    [ [ 1, -2, 1 ] ]
    gap> TriangulizedNullspaceMat(M);
    [ [ 1, -2, 1 ] ]
    gap> M:=[[1,2,3,1,1],[4,5,6,1,1],[7,8,9,1,1],[1,2,3,1,1]];
    [ [ 1, 2, 3, 1, 1 ], [ 4, 5, 6, 1, 1 ], [ 7, 8, 9, 1, 1 ], 
      [ 1, 2, 3, 1, 1 ] ]
    gap> NullspaceMat(M);
    [ [ 1, -2, 1, 0 ], [ -1, 0, 0, 1 ] ]
    gap> TriangulizedNullspaceMat(M);
    [ [ 1, 0, 0, -1 ], [ 0, 1, -1/2, -1/2 ] ]
    
    
    


  • characteristic polynomial:
    Please see section
    24.12.1 of the GAP reference manual
    for examples of characteristic polynomial of a
    square matrix (“CharacteristicPolynomial”) and
    section

    56.3
    for examples of the “characteristic polynomial”
    (called a “TracePolynomial”) of an
    element of a field extension.


  • character:
    GAP contains very extensive character theoretic functions
    and data libraries (including an interface the character table in the
    Atlas).
    Here is just one simple example.

    gap> G:=Group((1,2)(3,4),(1,2,3));
    Group([ (1,2)(3,4), (1,2,3) ])
    gap> T:=CharacterTable(G);
    CharacterTable( Alt( [ 1 .. 4 ] ) )
    gap> Display(T);
    CT1
    
         2  2  2  .  .
         3  1  .  1  1
    
           1a 2a 3a 3b
        2P 1a 1a 3b 3a
        3P 1a 2a 1a 1a
    
    X.1     1  1  1  1
    X.2     1  1  A /A
    X.3     1  1 /A  A
    X.4     3 -1  .  .
    
    A = E(3)^2
      = (-1-ER(-3))/2 = -1-b3
    gap> irr:=Irr(G);
    [ Character( CharacterTable( Alt( [ 1 .. 4 ] ) ), [ 1, 1, 1, 1 ] ),
      Character( CharacterTable( Alt( [ 1 .. 4 ] ) ), [ 1, 1, E(3)^2, E(3) ] ),
      Character( CharacterTable( Alt( [ 1 .. 4 ] ) ), [ 1, 1, E(3), E(3)^2 ] ),
      Character( CharacterTable( Alt( [ 1 .. 4 ] ) ), [ 3, -1, 0, 0 ] ) ]
    gap> Display(irr);
    [ [       1,       1,       1,       1 ],
      [       1,       1,  E(3)^2,    E(3) ],
      [       1,       1,    E(3),  E(3)^2 ],
      [       3,      -1,       0,       0 ] ]
    gap> chi:=irr[2]; gamma:=CG[3]; g:=Representative(gamma); g^chi;
    Character( CharacterTable( Alt( [ 1 .. 4 ] ) ), [ 1, 1, E(3)^2, E(3) ] )
    (1,2,3)^G
    (1,2,3)
    E(3)^2
    
    

    For further details and examples, see chapters
    69
    72 of the
    GAP reference manual.

  • brauer:
    Just a simple example of what GAP can do here.
    To construct a Brauer character table:

    gap> G:=Group((1,2)(3,4),(1,2,3));
    Group([ (1,2)(3,4), (1,2,3) ])
    gap> irr:=IrreducibleRepresentations(G,GF(7));
    [ [ (1,2)(3,4), (1,2,3) ] -> [ [ [ Z(7)^0 ] ], [ [ Z(7)^0 ] ] ],
      
    [ (1,2)(3,4), (1,2,3) ] -> [ [ [ Z(7)^0 ] ], [ [ Z(7)^4 ] ] ],
      
    [ (1,2)(3,4), (1,2,3) ] -> [ [ [ Z(7)^0 ] ], [ [ Z(7)^2 ] ] ],
      
    [ (1,2)(3,4), (1,2,3) ] -> [
          
    [ [ 0*Z(7), Z(7)^3, Z(7)^0 ], [ 0*Z(7), Z(7)^3, 0*Z(7) ], 
    [ Z(7)^0, Z(7)^3, 0*Z(7) ] ],
          [ [ 0*Z(7), Z(7)^0, 0*Z(7) ], 
    [ 0*Z(7), 0*Z(7), Z(7)^0 ], [ Z(7)^0, 0*Z(7), 0*Z(7) ] ]
         
    ] ]
    gap> brvals := List(irr,chi-> List(ConjugacyClasses(G),c->
    BrauerCharacterValue(Image(chi, Representative(c)))));
    [ [ 1, 1, 1, 1 ], [ 1, 1, E(3)^2, E(3) ], [ 1, 1, E(3), E(3)^2 ], 
    [ 3, -1, 0, 0 ] ]
    gap> Display(brvals);
    [ [       1,       1,       1,       1 ],
      
    [       1,       1,  E(3)^2,    E(3) ],
      
    [       1,       1,    E(3),  E(3)^2 ],
      
    [       3,      -1,       0,       0 ] ]
    gap>                                               
    

    List(ConjugacyClasses(G),c->BrauerCharacterValue(Image(chi, Representative(c)))));
    #Display(brvals);
    T:=CharacterTable(G);
    Display(T);
    –>


  • polynomial
    There are various ways to construct a polynomial in GAP.

    gap> Pts:=Z(7)^0*[1,2,3];
    [ Z(7)^0, Z(7)^2, Z(7) ]
    gap> Vals:=Z(7)^0*[1,2,6];
    [ Z(7)^0, Z(7)^2, Z(7)^3 ]
    gap> g:=InterpolatedPolynomial(GF(7),Pts,Vals);
    Z(7)^5*x_1^2+Z(7)
    

    Or:

    gap> p:=3;; F:=GF(p);;
    gap> R:=PolynomialRing(F,["x1","x2"]);
    PolynomialRing(..., [ x1, x2 ])
    gap> vars:=IndeterminatesOfPolynomialRing(R);;
    gap> x1:=vars[1]; x2:=vars[2];
    x1
    x2
    gap> p:=x1^5-x2^5;
    x1^5-x2^5
    gap> DivisorsMultivariatePolynomial(p,R);
    [ x1^4+x1^3*x2+x1^2*x2^2+x1*x2^3+x2^4, x1-x2 ]
    

    Or:

    gap> x:=X(Rationals);
    x_1
    gap> f:=x+x^2+1;
    x_1^2+x_1+1
    gap> Value(f,[x],[1]);
    3
    


  • factor
    To factor a polynomial in GAP, there is one command for
    univariate polynomials (“Factors”) and another command for
    multivariate polynomials (“DivisorsMultivariatePolynomial”).For a factoring a univariate polynomial,
    GAP provides only methods over finite fields
    and over subfields of cyclotomic fields.
    Please see the
    examples given in section

    64.10 “Polynomial Factorization”
    for more details.For multivariate polynomials,
    a very slow algorithm has been implemented in GAP
    and an interface to a very fast algorithm in
    Singular
    has been implemented for those who have both Singular and
    the GAP Singular package
    installed. The former of these was
    illustrated above in
    “polynomial” above.
    (Again, the ground field must be a finite field
    or a subfields of cyclotomic fields.)
    For the latter, please see the example
    in the (GAP-)Singular manual
    FactorsUsingSingularNC.


  • roots
    There are some situtations where GAP does find the roots
    of a polynomial but GAP does not do this generally.
    (The roots must generate either a finite field
    or a subfield of a cyclotomic field.) However, there is a package called

    RadiRoot
    which must be installed which does help to do this
    for polynomials with rational coefficients
    (radiroot itself requires other packages to be installed;
    please see the webpage for more details).The “Factors” command actually has an option which allows you to
    increase the groundfield so that a factorization actually
    returns the roots. Please see the
    examples given in section

    64.10 “Polynomial Factorization”
    for more details.Here is a second appoach.

    gap> p:=3; n:=4; F:=GF(p^n); c:=Random(F); r:=2;
    3
    4
    GF(3^4)
    Z(3^4)^79
    2
    gap>  x:=X(F,1); f:=x^r-c*x+c-1;
    x_1
    x_1^2+Z(3^4)^39*x_1+Z(3^4)^36
    gap>  F_f:=FieldExtension( F, f );
    AsField( GF(3^4), GF(3^8) )
    gap>  alpha:=RootOfDefiningPolynomial(F_f);
    Z(3^4)^36
    gap> Value(f,[x],[alpha]);
    0*Z(3)
    
    

    Here is a third. First, enter the following program:

    RootOfPolynomial:=function(f,R)
     local F0,Ff,a;
     F0:=CoefficientsRing(R);
     Ff:=FieldExtension(F0,f);
     a:=RootOfDefiningPolynomial(Ff);
     return a;
    end;
    

    Here’s how this can be used to find a root:

    gap> F:=Rationals;
    Rationals
    gap> x:=X(F,1); f:=x^2+x+1;
    x_1
    x_1^2+x_1+1
    gap> R:=PolynomialRing( F, [ x ]);
    PolynomialRing(..., [ x_1 ])
    gap> a:=RootOfPolynomial(f,R);
    E(3)
    gap> # check:
    gap> Value(f,[x],[a]);
    0
    

    Related links:

    1. In the GAP Forum:

      Hensel lifting discussion
      .
    2. In the manual,

      Galois groups
      .

  • evaluate:
    The relevant command is “Value”. There are several examples already on
    this page. For others, please see the examples given in section
    64.7 Multivariate polynomials of the manual.
    For sparse uivariate polynomials, there is also the command
    “ValuePol” in section
    23.6 of the manual.


  • integer power
    This is easy and intuitive:

    gap> a:=1000; n:=100000; m:=123;
    1000
    100000
    123
    gap> a^n mod m;
    1
    
    


  • matrix power:
    This too is easy and intuitive:

    gap> A:=[[1,2],[3,4]]; n:=100000; m:=123;
    [ [ 1, 2 ], [ 3, 4 ] ]
    100000
    123
    gap> A^n mod m;
    [ [ 1, 41 ], [ 0, 1 ] ]
    

  • polynomial power
    GAP allows you to do arithmetic over the polynomial
    ring R[x], where R = Z/nZ (where n is a positive integer).
    Here’s an example.

    gap> Z4:=ZmodnZ(4);
    (Integers mod 4)
    gap> R:=UnivariatePolynomialRing(Z4,1);
    PolynomialRing(..., [ x ])
    gap> x:=IndeterminatesOfPolynomialRing(R)[1];
    x
    gap> I:=TwoSidedIdealByGenerators( R,[x^8-x^0]);
    two-sided ideal in PolynomialRing(..., [ x ]), (1 generators)
    gap> gen:=x^8-x^0;
    x8-ZmodnZObj(1,4)
    gap> QuotientRemainder(R,x^8,gen);
    [ ZmodnZObj(1,4), ZmodnZObj(1,4) ]
    gap> QuotientRemainder(R,x^15,gen);
    [ x^7, x^7 ]
    gap> QuotientRemainder(R,x^15+x^8,gen);
    [ x^7+ZmodnZObj(1,4), x^7+ZmodnZObj(1,4) ]
    gap> PowerMod( R, x+x^0, 15, gen );
    ZmodnZObj(0,4)
    gap> PowerMod( R, x, 15, gen );
    x^7
    
    


  • Groebner basis
    GAP’s Groebner bases algorithms are relatively slow
    and are included mostly for simple examples and for
    teaching purposes. However, a GAP interface to a very
    fast algorithm in Singular
    has been implemented for those who have both Singular and
    the
    GAP Singular package
    installed. The former of these is
    illustrated in section
    64.17 Groebner bases of the GAP manual.
    For the latter, please see the example
    in the (GAP-)Singular manual
    GroebnerBasis.


  • normal subgroup:
    Here is an example:

    gap> G := AlternatingGroup( 5 );
    Group( (1,2,5), (2,3,5), (3,4,5) )
    gap> normal := NormalSubgroups( G );
    [ Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), [  ] ), 
      Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), 
        [ (1,2)(3,4), (1,3)(4,5), (1,4)(2,3) ] ) ]
    

    Related links:

    1. Please see Volkmar Felsch’s
      GAP Forum response to a related question.
    2. The

      xgap
      package displays subgroup lattices graphically.

  • abelian subgroup
    One idea to compute all the abelian subgroups is to compute all the
    subgroups then “filter” out the abelian ones.
    Here is an illustration, taked from a
    GAP Forum response Volkmar Felsch.

    gap> G := AlternatingGroup( 5 );
    Group( (1,2,5), (2,3,5), (3,4,5) )
    gap> classes := ConjugacyClassesSubgroups( G );
    [ ConjugacyClassSubgroups( Group( (1,2,5), (2,3,5), 
        (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), [  ] ) ),
      ConjugacyClassSubgroups( Group( (1,2,5), (2,3,5), 
        (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), 
        [ (2,3)(4,5) ] ) ), ConjugacyClassSubgroups( Group( (1,2,5), 
        (2,3,5), (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), 
        [ (3,4,5) ] ) ), ConjugacyClassSubgroups( Group( (1,2,5), 
        (2,3,5), (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), 
        [ (2,3)(4,5), (2,4)(3,5) ] ) ), ConjugacyClassSubgroups( Group( 
        (1,2,5), (2,3,5), (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), 
        (3,4,5) ), [ (1,2,3,4,5) ] ) ), ConjugacyClassSubgroups( Group( 
        (1,2,5), (2,3,5), (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), 
        (3,4,5) ), [ (3,4,5), (1,2)(4,5) ] ) ), 
      ConjugacyClassSubgroups( Group( (1,2,5), (2,3,5), 
        (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), 
        [ (1,2,3,4,5), (2,5)(3,4) ] ) ), ConjugacyClassSubgroups( Group( 
        (1,2,5), (2,3,5), (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), 
        (3,4,5) ), [ (2,3)(4,5), (2,4)(3,5), (3,4,5) ] ) ), 
      ConjugacyClassSubgroups( Group( (1,2,5), (2,3,5), (3,4,5) ), Group( 
        (1,2,5), (2,3,5), (3,4,5) ) ) ]
    gap> cl := classes[4];
    ConjugacyClassSubgroups( Group( (1,2,5), (2,3,5), 
    (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), 
    [ (2,3)(4,5), (2,4)(3,5) ] ) )
    gap> length := Size( cl );
    5
    gap> rep := Representative( cl );
    Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), 
    [ (2,3)(4,5), (2,4)(3,5) ] )
    gap> order := Size( rep );
    4
    gap> IsAbelian( rep );
    true
    gap> abel := Filtered( classes, cl -> IsAbelian( Representative( cl ) ) );
    [ ConjugacyClassSubgroups( Group( (1,2,5), (2,3,5), 
        (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), [  ] ) ),
      ConjugacyClassSubgroups( Group( (1,2,5), (2,3,5), 
        (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), 
        [ (2,3)(4,5) ] ) ), ConjugacyClassSubgroups( Group( (1,2,5), 
        (2,3,5), (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), 
        [ (3,4,5) ] ) ), ConjugacyClassSubgroups( Group( (1,2,5), 
        (2,3,5), (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), (3,4,5) ), 
        [ (2,3)(4,5), (2,4)(3,5) ] ) ), ConjugacyClassSubgroups( Group( 
        (1,2,5), (2,3,5), (3,4,5) ), Subgroup( Group( (1,2,5), (2,3,5), 
        (3,4,5) ), [ (1,2,3,4,5) ] ) ) ]
    

  • homology
    This depends on how the group is given. For example, suppose that
    G is a permutation group with generators genG and
    H is a permutation group with generators genH. To find a
    homomorphism from G to H, one may use the
    “GroupHomomorphismByImages” or “GroupHomomorphismByImagesNC”
    commands. For examples of the syntax, please see
    section
    38.1 Creating Group Homomorphisms.Here’s an illustration of how to convert a finitely presented
    group into a permutation group.

    gap> p:=7;
    7
    gap> G:=PSL(2,p);
    Group([ (3,7,5)(4,8,6), (1,2,6)(3,4,8) ])
    gap> H:=SchurCover(G);
    fp group of size 336 on the generators [ f1, f2, f3 ]
    gap> iso:=IsomorphismPermGroup(H);
    [ f1, f2, f3 ] -> [ (1,2,4,3)(5,9,7,10)(6,11,8,12)(13,14,15,16),
      (2,5,6)(3,7,8)(11,13,14)(12,15,16), (1,4)(2,3)(5,7)(6,8)(9,10)(11,12)(13,
        15)(14,16) ]
    gap> H0:=Image(iso);                       # 2-cover of PSL2
    Group([ (1,2,4,3)(5,9,7,10)(6,11,8,12)(13,14,15,16),
      (2,5,6)(3,7,8)(11,13,14)(12,15,16), (1,4)(2,3)(5,7)(6,8)(9,10)(11,12)(13,
        15)(14,16) ])
    gap> IdGroup(H0);
    [ 336, 114 ]
    gap> IdGroup(SL(2,7));
    [ 336, 114 ]
    gap>                
    

  • semi-direct product(Contributed by Nilo de Roock):
    As you can easily verify, D8 is isomorphic to C2:C4. Or in GAP…

    N:=CyclicGroup(IsPermGroup,4);
    G:=CyclicGroup(IsPermGroup,2);
    AutN:=AutomorphismGroup(N);
    f:=GroupHomomorphismByImages(G,AutN,GeneratorsOfGroup(G),[Elements(AutN)[2]]);
    NG:=SemidirectProduct(G,f,N);
    

    Verify with

    StructureDescription(NG);
    

  • semi-direct products(Contributed by Nilo de Roock):
    The following shows how to construct all non-abelian groups
    of order 12 as semi-direct products. These products are not
    trivial yet small enough to verify by hand.

    #D12 = (C2 x C2) : C3
    G1:=CyclicGroup(IsPermGroup,2);
    G2:=CyclicGroup(IsPermGroup,2);
    G:=DirectProduct(G1,G2);
    N:=CyclicGroup(IsPermGroup,3);
    AutN:=AutomorphismGroup(N);
    f:=GroupHomomorphismByImages(G,AutN,[Elements(G)[1],Elements(G)[2],Elements(G)[3],Elements(G)[4]],[Elements(AutN)[1],Elements(AutN)[2],Elements(AutN)[1],Elements(AutN)[2]]);
    NG:=SemidirectProduct(G,f,N);
    Print(str(NG));
    Print("\n");
    
    #T = C4 : C3
    G:=CyclicGroup(IsPermGroup,4);
    N:=CyclicGroup(IsPermGroup,3);
    AutN:=AutomorphismGroup(N);
    f:=GroupHomomorphismByImages(G,AutN,[Elements(G)[1],Elements(G)[2],Elements(G)[3],Elements(G)[4]],[Elements(AutN)[1],Elements(AutN)[2],Elements(AutN)[1],Elements(AutN)[2]]);
    NG:=SemidirectProduct(G,f,N);
    Print(str(NG));
    Print("\n");
    
    #A4 = C3 : (C2 x C2)
    G:=CyclicGroup(IsPermGroup,3);
    N1:=CyclicGroup(IsPermGroup,2);
    N2:=CyclicGroup(IsPermGroup,2);
    N:=DirectProduct(G1,G2);
    AutN:=AutomorphismGroup(N);
    f:=GroupHomomorphismByImages(G,AutN,[Elements(G)[1],Elements(G)[2],Elements(G)[3]],[Elements(AutN)[1],Elements(AutN)[4],Elements(AutN)[5]]);
    NG:=SemidirectProduct(G,f,N);
    Print(str(NG));
    Print("\n");
    

  • cohomology
    GAP will compute the Schur multiplier
    H2(G,C) using the
    “AbelianInvariantsMultiplier” command.
    Here is an example showing how to find H2(A5,C),
    where A5 is the alternating group on 5 letters.

    gap> A5:=AlternatingGroup(5);
    Alt( [ 1 .. 5 ] )
    gap> AbelianInvariantsMultiplier(A5);
    [ 2 ]
    

    So, H2(A5,C) is Z/2Z.

    Related links:

    1. See section

      37.23
      and
      section

      37.24
      of the GAP manual.
    2. See D. Holt’s GAP package
      cohomolo.

A water sharing problem and Sage

Here is the problem: There are n students, each with a cup that could contain water or could be empty. We assume that the cups have no limit and that there is at least one student with some water. Order the students from those who have the most water to those who have the least, with arbitrary choices in cases of ties. Following this order, each student takes turn sharing his water equally with the other students. In other words, a student with r ounces of water will give r/n ounces to each of the others (and hence have r/n ounces remaining for himself).

Q1: Is it possible that there is an initial configuration so that, at the end, all the students have exactly the same amount they started with?
This would mean that sharing was in fact not a loss.

Q2: Will this process “usually even out” the distribution of water?
(Part of the answer is to try to make this question more precise!)

Q3: Does the initial ordering of the students (ie, who pours when) matter?

This can be viewed as a matrix question. The action of the i-th student sharing water can be viewed as a linear transformation A_i on the vector of water quantities, indexed by the students. I think that in general the answer to Q1 is “yes.” Here is the Sage code in the case of 3 students:

sage: A1 = matrix(QQ, [[1/3,0,0],[1/3,1,0],[1/3,0,1]])
sage: A2 = matrix(QQ, [[1,1/3,0],[0,1/3,0],[0,1/3,1]])
sage: A3 = matrix(QQ, [[1,0,1/3],[0,1,1/3],[0,0,1/3]])
sage: (A1*A2*A3).eigenvectors_right()                 
[(1, [(1, 2, 3)], 1), 
 (0.1851851851851852? - 0.05237828008789241?*I, 
 [(1, 0.?e-57 + 1.414213562373095?*I, -1.000000000000000? - 1.414213562373095?*I)], 1), 
 (0.1851851851851852? + 0.05237828008789241?*I, 
 [(1, 0.?e-57 - 1.414213562373095?*I, -1.000000000000000? + 1.414213562373095?*I)], 1)]

In other words, the product A_1A_2A_3 has eigenvalue \lambda = 1 with eigenvector \vec{v} = (1,2,3). This means that in the case that student 1 has 3 ounces of water, student 2 has 2 ounces of water, student 3 has 1 ounce of water, if each student shares water as explained in the problem, at the end of the process, they will each have the same amount as they started with.

I don’t know the answer to Q2 or Q3.

SymPy and the GSoC

Google has announced the results for Google Summer of Code. The following projects have been accepted for SymPy:

(Project, Student, Mentor, Link to proposal on the wiki)
– Category Theory Module, Sergiu Ivanov, Tom Bachmann
– Density Operators for Quantum Module in sympy.physics.quantum, Guru
Devanla, Brian Granger (co-mentor Sean Vig)
– Enhancements to sympy.physics.mechanics, Angadh Nanjangud, Gilbert Gede
– Group Theory, Aleksandar Makelov, David Joyner (Aaron Meurer co-mentor)
– Implicit Plotting Module, Bharath M R, Aaron Meurer
– Vector Analysis, Stefan Krastanov, Matthew Rocklin

I will help mentor Aleksandar Makelov’s work on group theory. He is a freshman at Harvard.

Bifid cipher and Sage

The Bifid cipher was invented around 1901 by Felix Delastelle. It is a “fractional substitution” cipher, where letters are replaced by pairs of symbols from a smaller alphabet. The cipher uses a 5×5 square filled with some ordering of the alphabet, except that “i”‘s and “j”‘s are identified (this is a so-called Polybius square; there is a 6×6 analog if you add back in “j” and also append onto the usual 26 letter alphabet, the digits 0, 1, …, 9). According to Helen Gaines’ book “Cryptanalysis”, this type of cipher was used in the field by the German army during World War I.

The Bifid cipher was discusses in Alasdair McAndrew’s book on Cryptography and Sage. We shall follow his discussion. As an example of a Polybius square for the Bifid cipher, pick the key to be “encrypt” (as Alasdair does). In that case, the Polybius square is \left(\begin{array}{rrrrr}  E & N & C & R & Y \\  P & T & A & B & C \\  D & E & F & G & H \\  I & K & L & M & N \\  O & P & Q & R & S  \end{array}\right). BTW, the 6\times 6 analog is: \left(\begin{array}{rrrrrr}  E & N & C & R & Y & P \\  T & A & B & C & D & E \\  F & G & H & I & J & K \\  L & M & N & O & P & Q \\  R & S & T & U & V & W \\  X & Y & Z & 0 & 1 & 2  \end{array}\right).

Here is Sage code to produce the 6\times 6 case (the 5\times 5 case is in Alasdair’s book):

def bifid(pt, key):
    """
    INPUT:
        pt - plaintext string      (digits okay)
        key - short string for key (no repetitions, digits okay)
    
    OUTPUT:
        ciphertext from Bifid cipher (all caps, no spaces)

    This is the version of the Bifid cipher that uses the 6x6
    Polybius square.

    AUTHOR:
        Alasdair McAndrew

    EXAMPLES:
        sage: key = "encrypt"
        sage: pt = "meet me on monday at 8am"
        sage: bifid(pt, key)
        [[2, 5], [0, 0], [0, 0], [1, 0], [2, 5], [0, 0], [3, 0], 
         [0, 1], [2, 5], [3, 0], [0, 1], [1, 3], [1, 1], [0, 4], 
         [1, 1], [1, 0], [5, 4], [1, 1], [2, 5]]
        'HNHOKNTA5MEPEGNQZYG'

    """
    AS = AlphabeticStrings()
    A = list(AS.alphabet())+[str(x) for x in range(10)]
    # first make sure the letters are capitalized
    # and text has no spaces
    key0 = [x.capitalize() for x in key if not(x.isspace())]
    pt0 = [x.capitalize() for x in pt if not(x.isspace())]
    # create long key
    long_key = key0+[x for x in A if not(x in key0)]
    n = len(pt0)
    # the fractionalization
    pairs = [[long_key.index(x)//6, long_key.index(x)%6] for x in pt0]
    print pairs
    tmp_cipher = flatten([x[0] for x in pairs]+[x[1] for x in pairs])
    ct = join([long_key[6*tmp_cipher[2*i]+tmp_cipher[2*i+1]] for i in range(n)], sep="")
    return ct

def bifid_square(key):
    """
    Produced the Polybius square for the 6x6 Bifid cipher.

    EXAMPLE:
        sage: key = "encrypt"
        sage: bifid_square(key)
        [E N C R Y P]
        [T A B C D E]
        [F G H I J K]
        [L M N O P Q]
        [R S T U V W]
        [X Y Z 0 1 2]

    """
    AS = AlphabeticStrings()
    A = list(AS.alphabet())+[str(x) for x in range(10)]
    # first make sure the letters are capitalized
    # and text has no spaces
    key0 = [SR(x.capitalize()) for x in key if not(x.isspace())]
    # create long key
    long_key = key0+[SR(x) for x in A if not(x in key0)]
    # the fractionalization
    pairs = [[long_key.index(SR(x))//6, long_key.index(SR(x))%6] for x in A]
    f = lambda i,j: long_key[6*i+j] 
    M = matrix(SR, 6, 6, f)
    return M

Have fun!

uniform matroids and MDS codes

It is known that uniform (resp. paving) matroids correspond to MDS (resp. “almost MDS” codes). This post explains this connection.

An MDS code is an [n,k,d] linear error correcting block code C which meets the Singleton bound, d+k=n+1. A uniform matroid is a matroid for which all circuits are of size \geq r(M)+1, where r(M) is the rank of M. Recall, a circuit in a matroid M=(E,J) is a minimal dependent subset of E — that is, a dependent set whose proper subsets are all independent (i.e., all in J).

Consider a linear code C whose check matrix is an (n-k)\times n matrix H=(\vec{h}_1,\dots , \vec{h}_n). The vector matroid M=M[H] is a matroid for which the smallest sized dependency relation among the columns of H is determined by the check relations c_1\vec{h}_1 + \dots + c_n \vec{h}_n = H\vec{c}=\vec{0}, where \vec{c}=(c_1,\dots,c_n) is a codeword (in C which has minimum dimension d). Such a minimum dependency relation of H corresponds to a circuit of M=M[H].

Floyd-Warshall-Roy, 3

As in the previous post, let G=(V,E) be a weighted digraph having n vertices and let A=A_G denote itsn\times n adjacency matrix. We identify the vertices with the set \{1,2,\dots, n\}.

The previous post discussed the following result, due to Sturmfels et al.

Theorem: The entry of the matrix A\odot n-1 in row i and column j equals the length of a shortest path from vertex i to vertex j in G. (Here A\odot n-1 denotes the n-1-st tropical power of A.)

This post discusses an implementation in Python/Sage.

Consider the following class definition.

class TropicalNumbers:
    """
    Implements the tropical semiring.

    EXAMPLES:
        sage: T = TropicalNumbers()
        sage: print T
        Tropical Semiring

    """
    def __init__(self):
        self.identity = Infinity

    def __repr__(self):
        """
        Called to compute the "official" string representation of an object.
        If at all possible, this should look like a valid Python expression
        that could be used to recreate an object with the same value.

        EXAMPLES:
            sage: TropicalNumbers()
            TropicalNumbers()

        """
        return "TropicalNumbers()"

    def __str__(self):
        """
        Called to compute the "informal" string description of an object. 

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: print T
            Tropical Semiring

        """
        return "Tropical Semiring"

    def __call__(self, a):
        """
        Coerces a into the tropical semiring.

        EXAMPLES:
            sage: T(10)
            TropicalNumber(10)
            sage: print T(10)
            Tropical element 10 in Tropical Semiring
        """
        return TropicalNumber(a)

    def __contains__(self, a):
        """
        Implements "in".

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: a = T(10)
            sage: a in T
            True

        """
        if a in RR or a == Infinity:
            return a==Infinity or (RR(a) in RR)
        else:
            return a==Infinity or (RR(a.element) in RR)

class TropicalNumber:
    def __init__(self, a):
        self.element = a
        self.base_ring = TropicalNumbers()

    def __repr__(self):
        """
        Called to compute the "official" string representation of an object.
        If at all possible, this should look like a valid Python expression
        that could be used to recreate an object with the same value.

        EXAMPLES:

        """
        return "TropicalNumber(%s)"%self.element

    def __str__(self):
        """
        Called to compute the "informal" string description of an object. 

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: print T(10)
            Tropical element 10 in Tropical Semiring
        """
        return "%s"%(self.number())

    def number(self):
        return self.element

    def __add__(self, other):
        """
        Implements +. Assumes both self and other are instances of
        TropicalNumber class.

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: a = T(10)
            sage: a in T
            True
            sage: b = T(15)
            sage: a+b
            10

        """
        T = TropicalNumbers()
        return T(min(self.element,other.element))

    def __mul__(self, other):
        """
        Implements multiplication *.

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: a = T(10)
            sage: a in T
            True
            sage: b = T(15)
            sage: a*b
            25
        """
        T = TropicalNumbers()
        return T(self.element+other.element)

class TropicalMatrix:
    def __init__(self, A):
        T = TropicalNumbers()
        self.base_ring = T
        self.row_dimen = len(A)
        self.column_dimen = len(A[0])
        # now we coerce the entries into T
        A0 = A
        m = self.row_dimen
        n = self.column_dimen
        for i in range(m):
            for j in range(n):
                A0[i][j] = T(A[i][j])
        self.array = A0

    def matrix(self):
        """
        Returns the entries (as ordinary numbers).

        EXAMPLES:
            sage: A = [[0,1,3,7],[2,0,1,3],[4,5,0,1],[6,3,1,0]]
            sage: AT = TropicalMatrix(A)
            sage: AT.matrix()
            [[0, 1, 3, 7], [2, 0, 1, 3], [4, 5, 0, 1], [6, 3, 1, 0]]
        """
        m = self.row_dim()
        n = self.column_dim()
        A0 = [[0 for i in range(n)] for j in range(m)]
        for i in range(m):
            for j in range(n):
                A0[i][j] = (self.array[i][j]).number()
        return A0

    def row_dim(self):
        return self.row_dimen

    def column_dim(self):
        return self.column_dimen

    def __repr__(self):
        """
        Called to compute the "official" string representation of an object.
        If at all possible, this should look like a valid Python expression
        that could be used to recreate an object with the same value.

        EXAMPLES:

        """
        return "TropicalMatrix(%s)"%self.array

    def __str__(self):
        """
        Called to compute the "informal" string description of an object. 

        EXAMPLES:

        """
        return "Tropical matrix %s"%(self.matrix())

    def __add__(self, other):
        """
        Implements +. Assumes both self and other are instances of
        TropicalMatrix class.

        EXAMPLES:
            sage: A = [[1,2,Infinity],[3,Infinity,0]]
            sage: B = [[2,Infinity,1],[3,-1,1]]
            sage: AT = TropicalMatrix(A)
            sage: BT = TropicalMatrix(B)
            sage: AT
            TropicalMatrix([[TropicalNumber(1), TropicalNumber(2), TropicalNumber(+Infinity)],
              [TropicalNumber(3), TropicalNumber(+Infinity), TropicalNumber(0)]])
            sage: AT+BT
            [[TropicalNumber(1), TropicalNumber(2), TropicalNumber(1)],
             [TropicalNumber(3), TropicalNumber(-1), TropicalNumber(0)]]
        """
        A = self.array
        B = other.array
        C = []
        m = self.row_dim()
        n = self.column_dim()
        if m != other.row_dim:
            raise ValueError, "Row dimensions must be equal."
        if n != other.column_dim:
            raise ValueError, "Column dimensions must be equal."
        for i in range(m):
            row = [A[i][j]+B[i][j] for j in range(n)] # + as tropical numbers
            C.append(row)
        return C

    def __mul__(self, other):
        """
        Implements multiplication *.

        EXAMPLES:
            sage: A = [[1,2,Infinity],[3,Infinity,0]]
            sage: AT = TropicalMatrix(A)
            sage: B = [[2,Infinity],[-1,1],[Infinity,0]]
            sage: BT = TropicalMatrix(B)
            sage: AT*BT
             [[TropicalNumber(1), TropicalNumber(3)],
              [TropicalNumber(5), TropicalNumber(0)]]
            sage: A = [[0,1,3,7],[2,0,1,3],[4,5,0,1],[6,3,1,0]]
            sage: AT = TropicalMatrix(A)
            sage: A = [[0,1,3,7],[2,0,1,3],[4,5,0,1],[6,3,1,0]]
            sage: AT = TropicalMatrix(A)
            sage: print AT*AT*AT
            Tropical matrix [[0, 1, 2, 3], [2, 0, 1, 2], [4, 4, 0, 1], [5, 3, 1, 0]]
        """
        T = TropicalNumbers()
        A = self.matrix()
        B = other.matrix()
        C = []
        mA = self.row_dim()
        nA = self.column_dim()
        mB = other.row_dim()
        nB = other.column_dim()
        if nA != mB:
            raise ValueError, "Column dimension of A and row dimension of B must be equal."
        for i in range(mA):
            row = []
            for j in range(nB):
                c = T(Infinity)
                for k in range(nA):
                    c = c+T(A[i][k])*T(B[k][j])
                row.append(c.number())
            C.append(row)
        return TropicalMatrix(C)

This shows that the shortest distances of digraph with adjacency matrix \begin{pmatrix} 0&1&3&7\\ 2&0&1&3\\ 4&5&0&1\\ 6&3&1&0\end{pmatrix}, is equal to A\odot 3, which is equal to \begin{pmatrix} 0&1&2&3\\ 2&0&1&2\\ 4&4&0&1\\ 5&3&1&0\end{pmatrix}. This verifies an example given in chapter 1 of the book by Maclagan and Sturmfels, Introduction to Tropical Geometry .

Sage and the future of mathematics

I am not a biologist nor a chemist, and maybe neither are you, but suppose we were. Now, if I described a procedure, in “standard” detail, to produce a result XYZ, you would (based on your reasonably expected expertise in the field), follow the steps you believe were described and either reproduce XYZ or, if I was mistaken, not be able to reproduce XYZ. This is called scientific reproducibility. It is cructial to what I believe is one of the fundamental principles of science, namely Popper’s Falsifiability Criterion.

More and more people are arguing, correctly in my opinion, that in the computational realm, in particular in mathematical research which uses computational experiments, that much higher standards are needed. The Ars Technica article linked above suggests that “it’s time for a major revision of the scientific method.” Also, Victoria Stodden argues one must “facilitate reproducibility. Without this data may be open, but will remain de facto in the ivory tower.” The argument basically is that to reproduce computational mathematical experiments is unreasonably difficult, requiring more that a reasonable expertise. These days, it may in fact (unfortunately) require purchasing very expensive software, or possessing of very sophisticated programming skills, accessibility to special hardware, or (worse) guessing parameters and programming procedures only hinted at by the researcher.

Hopefully, Sage can play the role of a standard bearer for such computational reproducibility. Sage is free, open source and there is a publically available server it runs on (sagenb.org).

What government agencies should require such reproducibility? In my opinion, all scientific funding agencies (NSF, etc) should follow these higher standards of computational accountability.