--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
-- Copyright (C) 2007-2010, Gabriel Dos Reis.
-- All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

)abbrev domain PERMGRP PermutationGroup
++ Authors: G. Schneider, H. Gollan, J. Grabmeier
++ Date Created: 13 February 1987
++ Date Last Updated: 24 May 1991
++ Basic Operations:
++ Related Constructors: PermutationGroupExamples, Permutation
++ Also See: RepresentationTheoryPackage1
++ AMS Classifications:
++ Keywords: permutation, permutation group, group operation, word problem
++ References:
++   C. Sims: Determining the conjugacy classes of a permutation group,
++   in Computers in Algebra and Number Theory, SIAM-AMS Proc., Vol. 4,
++    Amer. Math. Soc., Providence, R. I., 1971, pp. 191-195
++ Description: 
++  PermutationGroup implements permutation groups acting
++  on a set S, i.e. all subgroups of the symmetric group of S,
++  represented as a list of permutations (generators). Note that
++  therefore the objects are not members of the \Language category
++  \spadtype{Group}.
++  Using the idea of base and strong generators by Sims,
++  basic routines and algorithms
++  are implemented so that the word problem for
++  permutation groups can be solved.
--++  Note: we plan to implement lattice operations on the subgroup
--++  lattice in a later release

PermutationGroup(S:SetCategory): public == private where

  L    ==> List
  PERM ==> Permutation
  FSET ==> Set
  I    ==> Integer
  NNI  ==> NonNegativeInteger
  V    ==> Vector
  B    ==> Boolean
  OUT   ==> OutputForm
  SYM  ==> Symbol
  REC  ==> Record ( orb : L NNI , svc : V I )
  REC2 ==> Record(order:NNI,sgset:L V NNI,_
             gpbase:L NNI,orbs:L REC,mp:L S,wd:L L NNI)
  REC3 ==> Record(elt:V NNI,lst:L NNI)
  REC4 ==> Record(bool:B,lst:L NNI)

  public ==> SetCategory with

    coerce           : %         -> L PERM S
      ++ coerce(gp) returns the generators of the group {\em gp}.
    generators           : %         -> L PERM S
      ++ generators(gp) returns the generators of the group {\em gp}.
    elt              : (%,NNI)   -> PERM S
      ++ elt(gp,i) returns the i-th generator of the group {\em gp}.
    random           : (%,I)     -> PERM S
      ++ random(gp,i) returns a random product of maximal i generators
      ++ of the group {\em gp}.
    random           : %         -> PERM S
      ++ random(gp) returns a random product of maximal 20 generators
      ++ of the group {\em gp}.
      ++ Note: {\em random(gp)=random(gp,20)}.
    order            : %         -> NNI
      ++ order(gp) returns the order of the group {\em gp}.
    degree           : %         -> NNI
      ++ degree(gp) returns the number of points moved by all permutations
      ++ of the group {\em gp}.
    base             : %         -> L S
      ++ base(gp) returns a base for the group {\em gp}.
    strongGenerators : %         -> L PERM S
      ++ strongGenerators(gp) returns strong generators for
      ++ the group {\em gp}.
    wordsForStrongGenerators      : %         -> L L NNI
      ++ wordsForStrongGenerators(gp) returns the words for the strong
      ++ generators of the group {\em gp} in the original generators of
      ++ {\em gp}, represented by their indices in the list, given by
      ++ {\em generators}.
    coerce           : L PERM S  -> %
      ++ coerce(ls) coerces a list of permutations {\em ls} to the group
      ++ generated by this list.
    permutationGroup          : L PERM S  -> %
      ++ permutationGroup(ls) coerces a list of permutations {\em ls} to
      ++ the group generated by this list.
    orbit            : (%,S)     -> FSET S
      ++ orbit(gp,el) returns the orbit of the element {\em el} under the
      ++ group {\em gp}, i.e. the set of all points gained by applying
      ++ each group element to {\em el}.
    orbits           : %         -> FSET FSET S
      ++ orbits(gp) returns the orbits of the group {\em gp}, i.e.
      ++ it partitions the (finite) of all moved points.
    orbit            : (%,FSET S)-> FSET FSET S
      ++ orbit(gp,els) returns the orbit of the unordered
      ++ set {\em els} under the group {\em gp}.
    orbit            : (%,L S)   -> FSET L S
      ++ orbit(gp,ls) returns the orbit of the ordered
      ++ list {\em ls} under the group {\em gp}.
      ++ Note: return type is L L S temporarily because FSET L S has an error.
      -- (GILT DAS NOCH?)
    member?          : (PERM S, %)-> B
      ++ member?(pp,gp) answers the question, whether the
      ++ permutation {\em pp} is in the group {\em gp} or not.
    wordInStrongGenerators : (PERM S, %)-> L NNI
      ++ wordInStrongGenerators(p,gp) returns the word for the
      ++ permutation p in the strong generators of the group {\em gp},
      ++ represented by the indices of the list, given by {\em strongGenerators}.
    wordInGenerators : (PERM S, %)-> L NNI
      ++ wordInGenerators(p,gp) returns the word for the permutation p
      ++ in the original generators of the group {\em gp},
      ++ represented by the indices of the list, given by {\em generators}.
    support      : %         -> FSET S
      ++ support(gp) returns the points moved by the group {\em gp}.
    <              : (%,%)     -> B
      ++ gp1 < gp2 returns true if and only if {\em gp1}
      ++ is a proper subgroup of {\em gp2}.
    <=             : (%,%)     -> B
      ++ gp1 <= gp2 returns true if and only if {\em gp1}
      ++ is a subgroup of {\em gp2}.
      ++ Note: because of a bug in the parser you have to call this
      ++ function explicitly by {\em gp1 <=$(PERMGRP S) gp2}.
      -- (GILT DAS NOCH?)
    initializeGroupForWordProblem : %   -> Void
      ++ initializeGroupForWordProblem(gp) initializes the group {\em gp}
      ++ for the word problem.
      ++ Notes: it calls the other function of this name with parameters
      ++ 0 and 1: {\em initializeGroupForWordProblem(gp,0,1)}.
      ++ Notes: (1) be careful: invoking this routine will destroy the
      ++ possibly information about your group (but will recompute it again)
      ++ (2) users need not call this function normally for the soultion of
      ++ the word problem.
    initializeGroupForWordProblem :(%,I,I) -> Void
      ++ initializeGroupForWordProblem(gp,m,n) initializes the group
      ++ {\em gp} for the word problem.
      ++ Notes: (1) with a small integer you get shorter words, but the
      ++ routine takes longer than the standard routine for longer words.
      ++ (2) be careful: invoking this routine will destroy the possibly stored
      ++ information about your group (but will recompute it again).
      ++ (3) users need not call this function normally for the soultion of
      ++ the word problem.

  private ==> add

    -- representation of the object:

    Rep  := Record ( gens : L PERM S , information : REC2 )

    -- import of domains and packages

    import Permutation S
    import OutputForm
    import Symbol
    import Void

  --first the local variables

    sgs               : L V NNI       := []
    baseOfGroup       : L NNI         := []
    sizeOfGroup       : NNI           := 1
    degree            : NNI           := 0
    gporb             : L REC         := []
    out               : L L V NNI     := []
    outword           : L L L NNI     := []
    wordlist          : L L NNI       := []
    basePoint         : NNI           := 0
    newBasePoint      : B             := true
    supp              : L S           := []
    ord               : NNI           := 1
    wordProblem       : B             := true

  --local functions first, signatures:

    shortenWord:(L NNI, %)->L NNI
    times:(V NNI, V NNI)->V NNI
    strip:(V NNI,REC,L V NNI,L L NNI)->REC3
    orbitInternal:(%,L S )->L L S
    inv: V NNI->V NNI
    ranelt:(L V NNI,L L NNI, I)->REC3
    testIdentity:V NNI->B
    pointList: %->L S
    orbitWithSvc:(L V NNI ,NNI )->REC
    cosetRep:(NNI ,REC ,L V NNI )->REC3
    bsgs1:(L V NNI,NNI,L L NNI,I,%,I)->NNI
    computeOrbits: I->L NNI
    reduceGenerators: I->Void
    bsgs:(%, I, I)->NNI
    initialize: %->FSET PERM S
    knownGroup?: %->Void
    subgroup:(%, %)->B
    memberInternal:(PERM S, %, B)->REC4

  --local functions first, implementations:

    shortenWord ( lw : L NNI , gp : % ) : L NNI ==
      -- tries to shorten a word in the generators by removing identities
      gpgens : L PERM S := coerce gp
      orderList : L NNI := [ order gen for gen in gpgens ]
      newlw : L NNI := copy lw
      for i in 1.. maxIndex orderList repeat
        if orderList.i = 1 then
          while member?(i,newlw) repeat
          -- removing the trivial element
            pos := position(i,newlw)
            newlw := delete(newlw,pos)
      flag : B := true
      while flag repeat
        actualLength : NNI := (maxIndex newlw) pretend NNI
        pointer := actualLength
        test := newlw.pointer
        anzahl : NNI := 1
        flag := false
        while pointer > 1 repeat
          pointer := ( pointer - 1 )::NNI
          if newlw.pointer ~= test then
            -- don't get a trivial element, try next
            test := newlw.pointer
            anzahl := 1
          else
            anzahl := anzahl + 1
            if anzahl = orderList.test then
              -- we have an identity, so remove it
              for i in (pointer+anzahl)..actualLength repeat
                newlw.(i-anzahl) := newlw.i
              newlw := first(newlw, (actualLength - anzahl) :: NNI)
              flag := true
              pointer := 1
      newlw

    times ( p : V NNI , q : V NNI ) : V NNI ==
      -- internal multiplication of permutations
      [ qelt(p,qelt(q,i)) for i in 1..degree ]

    strip(element:V NNI,orbit:REC,group:L V NNI,words:L L NNI) : REC3 ==
      -- strip an element into the stabilizer
      actelt         := element
      schreierVector := orbit.svc
      point          := orbit.orb.1
      outlist        := nil()$(L NNI)
      entryLessZero  : B := false
      while not entryLessZero repeat
        entry := schreierVector.(actelt.point)
        entryLessZero  := negative? entry
        if  not entryLessZero then
          actelt := times(group.entry, actelt)
          if wordProblem then outlist := append ( words.(entry::NNI) , outlist )
      [ actelt , reverse outlist ]

    orbitInternal ( gp : % , startList : L S ) : L L S ==
      orbitList : L L S := [ startList ]
      pos  : I := 1
      while not zero? pos  repeat
        gpset : L PERM S := gp.gens
        for gen in gpset repeat
          newList  := nil()$(L S)
          workList := orbitList.pos
          for j in #workList..1 by -1 repeat
            newList := cons (gen(workList.j) , newList )
          if not member?( newList , orbitList ) then
            orbitList := cons ( newList , orbitList )
            pos  := pos + 1
        pos := pos - 1
      reverse orbitList

    inv ( p : V NNI ) : V NNI ==
      -- internal inverse of a permutation
      q : V NNI := new(degree,0)$(V NNI)
      for i in 1..degree repeat q.(qelt(p,i)) := i
      q

    ranelt ( group : L V NNI , word : L L NNI , maxLoops : I ) : REC3 ==
      -- generate a "random" element
      numberOfGenerators    := # group
      randomInteger : I     := 1 + (random()$Integer rem numberOfGenerators)
      randomElement : V NNI := group.randomInteger
      words                 := nil()$(L NNI)
      if wordProblem then words := word.(randomInteger::NNI)
      if positive? maxLoops then
        numberOfLoops : I  := 1 + (random()$Integer rem maxLoops)
      else
        numberOfLoops : I := maxLoops
      while positive? numberOfLoops repeat
        randomInteger : I := 1 + (random()$Integer rem numberOfGenerators)
        randomElement := times ( group.randomInteger , randomElement )
        if wordProblem then words := append ( word.(randomInteger::NNI) , words)
        numberOfLoops := numberOfLoops - 1
      [ randomElement , words ]

    testIdentity ( p : V NNI ) : B ==
      -- internal test for identity
      for i in 1..degree repeat qelt(p,i) ~= i => return false
      true

    pointList(group : %) : L S ==
      s : FSET S :=  brace()   -- empty set !!
      for perm in group.gens repeat
        s := union(s, support perm)
      parts s

    orbitWithSvc ( group : L V NNI , point : NNI ) : REC ==
      -- compute orbit with Schreier vector, "-2" means not in the orbit,
      -- "-1" means starting point, the PI correspond to generators
      newGroup := nil()$(L V NNI)
      for el in group repeat
        newGroup := cons ( inv el , newGroup )
      newGroup               := reverse newGroup
      orbit          : L NNI := [ point ]
      schreierVector : V I   := new ( degree , -2 )
      schreierVector.point   := -1
      position : I := 1
      while not zero? position repeat
        for i in 1..#newGroup repeat
          newPoint := orbit.position
          newPoint := newGroup.i.newPoint
          if not member? ( newPoint , orbit ) then
            orbit                   := cons ( newPoint , orbit )
            position                := position + 1
            schreierVector.newPoint := i
        position := position - 1
      [ reverse orbit , schreierVector ]

    cosetRep ( point : NNI , o : REC , group : L V NNI ) : REC3 ==
      ppt          := point
      xelt : V NNI := [ n for n in 1..degree ]
      word         := nil()$(L NNI)
      oorb         := o.orb
      osvc         := o.svc
      while positive? degree repeat
        p := osvc.ppt
        negative? p => return [ xelt , word ]
        x    := group.p
        xelt := times ( x , xelt )
        if wordProblem then word := append ( wordlist.p , word )
        ppt  := x.ppt
      [xelt,word]

    bsgs1 (group:L V NNI,number1:NNI,words:L L NNI,maxLoops:I,gp:%,diff:I)_
        : NNI ==
      -- try to get a good approximation for the strong generators and base
      ort: REC
      k1: NNI
      i : NNI
      for i: free in number1..degree repeat
        ort := orbitWithSvc ( group , i )
        k   := ort.orb
        k1  := # k
        if not one? k1 then leave
      gpsgs := nil()$(L V NNI)
      words2 := nil()$(L L NNI)
      gplength : NNI := #group
      jj: NNI
      for jj: free in 1..gplength repeat if (group.jj).i ~= i then leave
      for k in 1..gplength repeat
        el2 := group.k
        if el2.i ~= i then
          gpsgs := cons ( el2 , gpsgs )
          if wordProblem then words2 := cons ( words.k , words2 )
        else
          gpsgs := cons ( times ( group.jj , el2 ) , gpsgs )
          if wordProblem _
            then words2 := cons ( append ( words.jj , words.k ) , words2 )
      group2 := nil()$(L V NNI)
      words3 := nil()$(L L NNI)
      j : I  := 15
      while positive? j repeat
        -- find generators for the stabilizer
        ran := ranelt ( group , words , maxLoops )
        str := strip ( ran.elt , ort , group , words )
        el2 := str.elt
        if not testIdentity el2 then
          if not member?(el2,group2) then
            group2 := cons ( el2 , group2 )
            if wordProblem then
              help : L NNI := append ( reverse str.lst , ran.lst )
              help         := shortenWord ( help , gp )
              words3       := cons ( help , words3 )
            j := j - 2
        j := j - 1
      -- this is for word length control
      if wordProblem then maxLoops    := maxLoops - diff
      if ( null group2 ) or negative? maxLoops then
        sizeOfGroup := k1
        baseOfGroup := [ i ]
        out         := [ gpsgs ]
        outword     := [ words2 ]
        return sizeOfGroup
      k2          := bsgs1 ( group2 , i + 1 , words3 , maxLoops , gp , diff )
      sizeOfGroup := k1 * k2
      out         := append ( out , [ gpsgs ] )
      outword     := append ( outword , [ words2 ] )
      baseOfGroup := cons ( i , baseOfGroup )
      sizeOfGroup

    computeOrbits ( kkk : I ) : L NNI ==
      -- compute the orbits for the stabilizers
      sgs         := nil()
      orbitLength := nil()$(L NNI)
      gporb       := nil()
      for i in 1..#baseOfGroup repeat
        sgs         := append ( sgs , out.i )
        pt          := #baseOfGroup - i + 1
        obs         := orbitWithSvc ( sgs , baseOfGroup.pt )
        orbitLength := cons ( #obs.orb , orbitLength )
        gporb       := cons ( obs , gporb )
      gporb := reverse gporb
      reverse orbitLength

    reduceGenerators ( kkk : I ) : Void ==
      -- try to reduce number of strong generators
      orbitLength := computeOrbits ( kkk )
      sgs         := nil()
      wordlist    := nil()
      for i in 1..(kkk-1) repeat
        sgs := append ( sgs , out.i )
        if wordProblem then wordlist := append ( wordlist , outword.i )
      removedGenerator := false
      baseLength : NNI := #baseOfGroup
      for nnn in kkk..(baseLength-1) repeat
        sgs := append ( sgs , out.nnn )
        if wordProblem then wordlist := append ( wordlist , outword.nnn )
        pt  := baseLength - nnn + 1
        obs := orbitWithSvc ( sgs , baseOfGroup.pt )
        i : NNI  := 1
        while not ( i > # out.nnn ) repeat
          pos  := position ( out.nnn.i , sgs )
          sgs2 := delete(sgs, pos)
          obs2 := orbitWithSvc ( sgs2 , baseOfGroup.pt )
          if # obs2.orb = orbitLength.nnn then
            test := true
            for j in (nnn+1)..(baseLength-1) repeat
              pt2  := baseLength - j + 1
              sgs2 := append ( sgs2 , out.j )
              obs2 := orbitWithSvc ( sgs2 , baseOfGroup.pt2 )
              if # obs2.orb ~= orbitLength.j then
                test := false
                leave
            if test then
              removedGenerator := true
              sgs              := delete (sgs, pos)
              if wordProblem then wordlist    := delete(wordlist, pos)
              out.nnn          := delete (out.nnn, i)
              if wordProblem then _
                outword.nnn := delete(outword.nnn, i )
            else
              i := i + 1
          else
            i := i + 1
      if removedGenerator then orbitLength := computeOrbits ( kkk )


    bsgs ( group : % ,maxLoops : I , diff : I ) : NNI ==
      -- the MOST IMPORTANT part of the package
      supp   := pointList group
      degree := # supp
      if degree = 0 then
        sizeOfGroup := 1
        sgs         := [ [ 0 ] ]
        baseOfGroup := nil()
        gporb       := nil()
        return sizeOfGroup
      newGroup := nil()$(L V NNI)
      gp       : L PERM S := group.gens
      words := nil()$(L L NNI)
      for ggg in 1..#gp repeat
        q := new(degree,0)$(V NNI)
        for i in 1..degree repeat
          newEl := elt(gp.ggg,supp.i)
          pos2  := position ( newEl , supp )
          q.i   := pos2 pretend NNI
        newGroup := cons ( q , newGroup )
        if wordProblem then words    := cons(list ggg, words)
      if maxLoops < 1 then
        -- try to get the (approximate) base length
        if zero? (# ((group.information).gpbase)) then
          wordProblem := false
          k           := bsgs1 ( newGroup , 1 , words , 20 , group , 0 )
          wordProblem := true
          maxLoops    := (# baseOfGroup) - 1
        else
          maxLoops    := (# ((group.information).gpbase)) - 1
      k       := bsgs1 ( newGroup , 1 , words , maxLoops , group , diff )
      kkk : I := 1
      newGroup := reverse newGroup
      noAnswer : B := true
      z: V NNI
      while noAnswer repeat
        reduceGenerators kkk
-- *** Here is former "bsgs2" *** --
        -- test whether we have a base and a strong generating set
        sgs := nil()
        wordlist := nil()
        for i in 1..(kkk-1) repeat
          sgs := append ( sgs , out.i )
          if wordProblem then wordlist := append ( wordlist , outword.i )
        noresult : B := true
        word3: L NNI
        word: L NNI
        for i in kkk..#baseOfGroup while noresult repeat
          sgs    := append ( sgs , out.i )
          if wordProblem then wordlist := append ( wordlist , outword.i )
          gporbi := gporb.i
          for pt in gporbi.orb while noresult repeat
            ppp   := cosetRep ( pt , gporbi , sgs )
            y1    := inv ppp.elt
            word3 := ppp.lst
            for jjj in 1..#sgs while noresult repeat
              word         := nil()$(L NNI)
              z            := times ( sgs.jjj , y1 )
              if wordProblem then word := append ( wordlist.jjj , word )
              ppp          := cosetRep ( (sgs.jjj).pt , gporbi , sgs )
              z            := times ( ppp.elt , z )
              if wordProblem then word := append ( ppp.lst , word )
              newBasePoint := false
              for j in (i-1)..1 by -1 while noresult repeat
                s := gporb.j.svc
                p := gporb.j.orb.1
                while positive? degree and noresult repeat
                  entry := s.(z.p)
                  if negative? entry then
                    if entry = -1 then leave
                    basePoint := j::NNI
                    noresult := false
                  else
                    ee := sgs.entry
                    z  := times ( ee , z )
                    if wordProblem then word := append ( wordlist.entry , word )
              if noresult then
                basePoint    := 1
                newBasePoint := true
                noresult := testIdentity z
        noAnswer := not (testIdentity z)
        if noAnswer then
          -- we have missed something
          word2 := nil()$(L NNI)
          if wordProblem then
            for wd in word3 repeat
              ttt := newGroup.wd
              while not (testIdentity ttt) repeat
                word2 := cons ( wd , word2 )
                ttt   := times ( ttt , newGroup.wd )
            word := append ( word , word2 )
            word := shortenWord ( word , group )
          if newBasePoint then
            for i in 1..degree repeat
              if z.i ~= i then
                baseOfGroup := append ( baseOfGroup , [ i ] )
                leave
            out := cons (list  z, out )
            if wordProblem then outword := cons (list word , outword )
          else
            out.basePoint := cons ( z , out.basePoint )
            if wordProblem then outword.basePoint := cons(word ,outword.basePoint )
          kkk := basePoint
      sizeOfGroup  := 1
      for j in 1..#baseOfGroup repeat
        sizeOfGroup := sizeOfGroup * # gporb.j.orb
      sizeOfGroup


    initialize ( group : % ) : FSET PERM S ==
      group2 := brace()$(FSET PERM S)
      gp : L PERM S := group.gens
      for gen in gp repeat
        if positive? degree gen then insert!(gen, group2)
      group2

    knownGroup? (gp : %) : Void ==
      -- do we know the group already?
      result := gp.information
      if result.order = 0 then
        wordProblem       := false
        ord               := bsgs ( gp , 20 , 0 )
        result            := [ ord , sgs , baseOfGroup , gporb , supp , [] ]
        gp.information    := result
      else
        ord         := result.order
        sgs         := result.sgset
        baseOfGroup := result.gpbase
        gporb       := result.orbs
        supp        := result.mp
        wordlist    := result.wd

    subgroup ( gp1 : % , gp2 : % ) : B ==
      gpset1 := initialize gp1
      gpset2 := initialize gp2
      empty? difference (gpset1, gpset2) => true
      for el in parts gpset1 repeat
        not member? (el,  gp2) => return false
      true

    memberInternal ( p : PERM S , gp : % , flag : B ) : REC4 ==
      -- internal membership testing
      supp     := pointList gp
      outlist := nil()$(L NNI)
      mP : L S := parts support p
      for x in mP repeat
        not member? (x, supp) => return [ false , nil()$(L NNI) ]
      if flag then
        member? ( p , gp.gens ) => return [ true , nil()$(L NNI) ]
        knownGroup? gp
      else
        result := gp.information
        if #(result.wd) = 0 then
          initializeGroupForWordProblem gp
        else
          ord         := result.order
          sgs         := result.sgset
          baseOfGroup := result.gpbase
          gporb       := result.orbs
          supp        := result.mp
          wordlist    := result.wd
      degree := # supp
      pp := new(degree,0)$(V NNI)
      for i in 1..degree repeat
        el   := p(supp.i)
        pos  := position ( el , supp )
        pp.i := pos::NNI
      words := nil()$(L L NNI)
      if wordProblem then
        for i in 1..#sgs repeat
          lw : L NNI := [ (#sgs - i + 1)::NNI ]
          words := cons ( lw , words )
      for i in #baseOfGroup..1 by -1 repeat
        str := strip ( pp , gporb.i , sgs , words )
        pp := str.elt
        if wordProblem then outlist := append ( outlist , str.lst )
      [ testIdentity pp , reverse outlist ]

  --now the exported functions

    coerce ( gp : % ) : L PERM S == gp.gens
    generators ( gp : % ) : L PERM S == gp.gens

    strongGenerators ( group ) ==
      knownGroup? group
      degree := # supp
      strongGens := nil()$(L PERM S)
      for i in sgs repeat
        pairs := nil()$(L L S)
        for j in 1..degree repeat
          pairs := cons ( [ supp.j , supp.(i.j) ] , pairs )
        strongGens := cons ( coerceListOfPairs pairs , strongGens )
      reverse strongGens

    elt ( gp , i ) == (gp.gens).i

    support ( gp ) == brace pointList gp

    random ( group , maximalNumberOfFactors ) ==
      maximalNumberOfFactors < 1 => 1$(PERM S)
      gp : L PERM S := group.gens
      numberOfGenerators := # gp
      randomInteger : I  := 1 + (random()$Integer rem numberOfGenerators)
      randomElement      := gp.randomInteger
      numberOfLoops : I  := 1 + (random()$Integer rem maximalNumberOfFactors)
      while positive? numberOfLoops repeat
        randomInteger : I  := 1 + (random()$Integer rem numberOfGenerators)
        randomElement := gp.randomInteger * randomElement
        numberOfLoops := numberOfLoops - 1
      randomElement

    random ( group ) == random ( group , 20 )

    order ( group ) ==
      knownGroup? group
      ord

    degree ( group ) == # pointList group

    base ( group ) ==
      knownGroup? group
      groupBase := nil()$(L S)
      for i in baseOfGroup repeat
        groupBase := cons ( supp.i , groupBase )
      reverse groupBase

    wordsForStrongGenerators ( group ) ==
      knownGroup? group
      wordlist

    coerce ( gp : L PERM S ) : % ==
      result : REC2 := [ 0 , [] , [] , [] , [] , [] ]
      group         := [ gp , result ]

    permutationGroup ( gp : L PERM S ) : % ==
      result : REC2 := [ 0 , [] , [] , [] , [] , [] ]
      group         := [ gp , result ]

    coerce(group: %) : OUT ==
      outList := nil()$(L OUT)
      gp : L PERM S := group.gens
      for i in (maxIndex gp)..1 by -1 repeat
        outList := cons(coerce gp.i, outList)
      postfix(outputForm(">":SYM),postfix(commaSeparate outList,outputForm("<":SYM)))

    orbit ( gp : % , el : S ) : FSET S ==
      elList : L S := [ el ]
      outList      := orbitInternal ( gp , elList )
      outSet       := brace()$(FSET S)
      for i in 1..#outList repeat
        insert! ( outList.i.1 , outSet )
      outSet

    orbits ( gp ) ==
      spp    := support gp
      orbits := nil()$(L FSET S)
      while positive? cardinality spp repeat
        el       := extract! spp
        orbitSet := orbit ( gp , el )
        orbits   := cons ( orbitSet , orbits )
        spp      := difference ( spp , orbitSet )
      brace orbits

    member? (p,  gp) ==
      wordProblem := false
      mi := memberInternal ( p , gp , true )
      mi.bool

    wordInStrongGenerators (p, gp ) ==
      mi := memberInternal ( inv p , gp , false )
      not mi.bool => error "p is not an element of gp"
      mi.lst

    wordInGenerators (p,  gp) ==
      lll : L NNI := wordInStrongGenerators (p, gp)
      outlist := nil()$(L NNI)
      for wd in lll repeat
        outlist := append ( outlist , wordlist.wd )
      shortenWord ( outlist , gp )

    gp1 < gp2 ==
      not empty? difference ( support gp1 , support gp2 ) => false
      not subgroup ( gp1 , gp2 ) => false
      order gp1 = order gp2 => false
      true

    gp1 <= gp2 ==
      not empty?  difference ( support gp1 , support gp2 ) => false
      subgroup ( gp1 , gp2 )

    gp1 = gp2 ==
      support gp1 ~= support gp2 => false
      if #(gp1.gens) <= #(gp2.gens) then
        not subgroup ( gp1 , gp2 ) => return false
      else
        not subgroup ( gp2 , gp1 ) => return false
      order gp1 = order gp2 => true
      false

    orbit ( gp : % , startSet : FSET S ) : FSET FSET S ==
      startList : L S := parts startSet
      outList         := orbitInternal ( gp , startList )
      outSet          := brace()$(FSET FSET S)
      for i in 1..#outList repeat
        newSet : FSET S := brace outList.i
        insert! ( newSet , outSet )
      outSet

    orbit ( gp : % , startList : L S ) : FSET L S ==
      brace orbitInternal(gp, startList)

    initializeGroupForWordProblem ( gp , maxLoops , diff ) ==
      wordProblem    := true
      ord            := bsgs ( gp , maxLoops , diff )
      gp.information := [ ord , sgs , baseOfGroup , gporb , supp , wordlist ]

    initializeGroupForWordProblem ( gp ) == initializeGroupForWordProblem ( gp , 0 , 1 )

)abbrev package PGE PermutationGroupExamples
++ Authors: M. Weller, G. Schneider, J. Grabmeier
++ Date Created: 20 February 1990
++ Date Last Updated: 09 June 1990
++ Basic Operations:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++  J. Conway, R. Curtis, S. Norton, R. Parker, R. Wilson:
++   Atlas of Finite Groups, Oxford, Clarendon Press, 1987
++ Description: 
++   PermutationGroupExamples provides permutation groups for
++   some classes of groups: symmetric, alternating, dihedral, cyclic,
++   direct products of cyclic, which are in fact the finite abelian groups
++   of symmetric groups called Young subgroups.
++   Furthermore, Rubik's group as permutation group of 48 integers and a list
++   of sporadic simple groups derived from the atlas of finite groups.

PermutationGroupExamples():public == private where

    L          ==> List
    I          ==> Integer
    PI         ==> PositiveInteger
    NNI        ==> NonNegativeInteger
    PERM       ==> Permutation
    PERMGRP   ==> PermutationGroup

    public ==> with

      symmetricGroup:       PI        -> PERMGRP I
        ++ symmetricGroup(n) constructs the symmetric group {\em Sn}
        ++ acting on the integers 1,...,n, generators are the
        ++ {\em n}-cycle {\em (1,...,n)} and the 2-cycle {\em (1,2)}.
      symmetricGroup:       L I       -> PERMGRP I
        ++ symmetricGroup(li) constructs the symmetric group acting on
        ++ the integers in the list {\em li}, generators are the
        ++ cycle given by {\em li} and the 2-cycle {\em (li.1,li.2)}.
        ++ Note: duplicates in the list will be removed.
      alternatingGroup:     PI        -> PERMGRP I
        ++ alternatingGroup(n) constructs the alternating group {\em An}
        ++ acting on the integers 1,...,n,  generators are in general the
        ++ {\em n-2}-cycle {\em (3,...,n)} and the 3-cycle {\em (1,2,3)}
        ++ if n is odd and the product of the 2-cycle {\em (1,2)} with
        ++ {\em n-2}-cycle {\em (3,...,n)} and the 3-cycle {\em (1,2,3)}
        ++ if n is even.
      alternatingGroup:     L I       -> PERMGRP I
        ++ alternatingGroup(li) constructs the alternating group acting
        ++ on the integers in the list {\em li}, generators are in general the
        ++ {\em n-2}-cycle {\em (li.3,...,li.n)} and the 3-cycle
        ++ {\em (li.1,li.2,li.3)}, if n is odd and
        ++ product of the 2-cycle {\em (li.1,li.2)} with
        ++ {\em n-2}-cycle {\em (li.3,...,li.n)} and the 3-cycle
        ++ {\em (li.1,li.2,li.3)}, if n is even.
        ++ Note: duplicates in the list will be removed.
      abelianGroup:         L PI      -> PERMGRP I
        ++ abelianGroup([n1,...,nk]) constructs the abelian group that
        ++ is the direct product of cyclic groups with order {\em ni}.
      cyclicGroup:          PI        -> PERMGRP I
        ++ cyclicGroup(n) constructs the cyclic group of order n acting
        ++ on the integers 1,...,n.
      cyclicGroup:          L I       -> PERMGRP I
        ++ cyclicGroup([i1,...,ik]) constructs the cyclic group of
        ++ order k acting on the integers {\em i1},...,{\em ik}.
        ++ Note: duplicates in the list will be removed.
      dihedralGroup:        PI        -> PERMGRP I
        ++ dihedralGroup(n) constructs the dihedral group of order 2n
        ++ acting on integers 1,...,N.
      dihedralGroup:        L I       -> PERMGRP I
        ++ dihedralGroup([i1,...,ik]) constructs the dihedral group of
        ++ order 2k acting on the integers out of {\em i1},...,{\em ik}.
        ++ Note: duplicates in the list will be removed.
      mathieu11:            L I       -> PERMGRP I
        ++ mathieu11(li) constructs the mathieu group acting on the 11
        ++ integers given in the list {\em li}.
        ++ Note: duplicates in the list will be removed.
        ++ error, if {\em li} has less or more than 11 different entries.
      mathieu11:            ()        -> PERMGRP I
        ++ mathieu11 constructs the mathieu group acting on the
        ++ integers 1,...,11.
      mathieu12:            L I       -> PERMGRP I
        ++ mathieu12(li) constructs the mathieu group acting on the 12
        ++ integers given in the list {\em li}.
        ++ Note: duplicates in the list will be removed
        ++ Error: if {\em li} has less or more than 12 different entries.
      mathieu12:            ()        -> PERMGRP I
        ++ mathieu12 constructs the mathieu group acting on the
        ++ integers 1,...,12.
      mathieu22:            L I       -> PERMGRP I
        ++ mathieu22(li) constructs the mathieu group acting on the 22
        ++ integers given in the list {\em li}.
        ++ Note: duplicates in the list will be removed.
        ++ Error: if {\em li} has less or more than 22 different entries.
      mathieu22:            ()        -> PERMGRP I
        ++ mathieu22 constructs the mathieu group acting on the
        ++ integers 1,...,22.
      mathieu23:            L I       -> PERMGRP I
        ++ mathieu23(li) constructs the mathieu group acting on the 23
        ++ integers given in the list {\em li}.
        ++ Note: duplicates in the list will be removed.
        ++ Error: if {\em li} has less or more than 23 different entries.
      mathieu23:            ()        -> PERMGRP I
        ++ mathieu23 constructs the mathieu group acting on the
        ++ integers 1,...,23.
      mathieu24:            L I       -> PERMGRP I
        ++ mathieu24(li) constructs the mathieu group acting on the 24
        ++ integers given in the list {\em li}.
        ++ Note: duplicates in the list will be removed.
        ++ Error: if {\em li} has less or more than 24 different entries.
      mathieu24:            ()        -> PERMGRP I
        ++ mathieu24 constructs the mathieu group acting on the
        ++ integers 1,...,24.
      janko2:               L I       -> PERMGRP I
        ++ janko2(li) constructs the janko group acting on the 100
        ++ integers given in the list {\em li}.
        ++ Note: duplicates in the list will be removed.
        ++ Error: if {\em li} has less or more than 100 different entries
      janko2:               ()        -> PERMGRP I
        ++ janko2 constructs the janko group acting on the
        ++ integers 1,...,100.
      rubiksGroup:          ()        -> PERMGRP I
        ++ rubiksGroup constructs the permutation group representing
        ++ Rubic's Cube acting on integers {\em 10*i+j} for
        ++ {\em 1 <= i <= 6}, {\em 1 <= j <= 8}.
        ++ The faces of Rubik's Cube are labelled in the obvious way
        ++ Front, Right, Up, Down, Left, Back and numbered from 1 to 6
        ++ in this given ordering, the pieces on each face
        ++ (except the unmoveable center piece) are clockwise numbered
        ++ from 1 to 8 starting with the piece in the upper left
        ++ corner. The moves of the cube are represented as permutations
        ++ on these pieces, represented as a two digit
        ++ integer {\em ij} where i is the numer of theface (1 to 6)
        ++ and j is the number of the piece on this face.
        ++ The remaining ambiguities are resolved by looking
        ++ at the 6 generators, which represent a 90 degree turns of the
        ++ faces, or from the following pictorial description.
        ++ Permutation group representing Rubic's Cube acting on integers
        ++ 10*i+j for 1 <= i <= 6, 1 <= j <=8.
        ++
        ++ \begin{verbatim}
        ++ Rubik's Cube:   +-----+ +-- B   where: marks Side # :
        ++                / U   /|/
        ++               /     / |         F(ront)    <->    1
        ++       L -->  +-----+ R|         R(ight)    <->    2
        ++              |     |  +         U(p)       <->    3
        ++              |  F  | /          D(own)     <->    4
        ++              |     |/           L(eft)     <->    5
        ++              +-----+            B(ack)     <->    6
        ++                 ^
        ++                 |
        ++                 D
        ++
        ++ The Cube's surface:
        ++                                The pieces on each side
        ++             +---+              (except the unmoveable center
        ++             |567|              piece) are clockwise numbered
        ++             |4U8|              from 1 to 8 starting with the
        ++             |321|              piece in the upper left
        ++         +---+---+---+          corner (see figure on the
        ++         |781|123|345|          left).  The moves of the cube
        ++         |6L2|8F4|2R6|          are represented as
        ++         |543|765|187|          permutations on these pieces.
        ++         +---+---+---+          Each of the pieces is
        ++             |123|              represented as a two digit
        ++             |8D4|              integer ij where i is the
        ++             |765|              # of the side ( 1 to 6 for
        ++             +---+              F to B (see table above ))
        ++             |567|              and j is the # of the piece.
        ++             |4B8|
        ++             |321|
        ++             +---+
        ++ \end{verbatim}
      youngGroup:           L I      -> PERMGRP I
        ++ youngGroup([n1,...,nk]) constructs the direct product of the
        ++ symmetric groups {\em Sn1},...,{\em Snk}.
      youngGroup:    Partition        -> PERMGRP I
        ++ youngGroup(lambda) constructs the direct product of the symmetric
        ++ groups given by the parts of the partition {\em lambda}.

    private ==> add

      -- import the permutation and permutation group domains:

      import PERM I
      import PERMGRP I

      -- import the needed map function:

      import ListFunctions2(L L I,PERM I)
      -- the internal functions:

      llli2gp(l:L L L I):PERMGRP I ==
        --++ Converts an list of permutations each represented by a list
        --++ of cycles ( each of them represented as a list of Integers )
        --++ to the permutation group generated by these permutations.
        (map(cycles,l))::PERMGRP I

      li1n(n:I):L I ==
        --++ constructs the list of integers from 1 to n
        [i for i in 1..n]

      -- definition of the exported functions:
      youngGroup(l:L I):PERMGRP I ==
        gens:= nil()$(L L L I)
        element:I:= 1
        for n in l | n > 1 repeat
          gens:=cons(list [i for i in element..(element+n-1)], gens)
          if n >= 3 then gens := cons([[element,element+1]],gens)
          element:=element+n
        llli2gp
          #gens = 0 => [[[1]]]
          gens

      youngGroup(lambda : Partition):PERMGRP I ==
        youngGroup(lambda::L(PI) pretend L I)

      rubiksGroup():PERMGRP I ==
        -- each generator represents a 90 degree turn of the appropriate
        -- side.
        f:L L I:=
         [[11,13,15,17],[12,14,16,18],[51,31,21,41],[53,33,23,43],[52,32,22,42]]
        r:L L I:=
         [[21,23,25,27],[22,24,26,28],[13,37,67,43],[15,31,61,45],[14,38,68,44]]
        u:L L I:=
         [[31,33,35,37],[32,34,36,38],[13,51,63,25],[11,57,61,23],[12,58,62,24]]
        d:L L I:=
         [[41,43,45,47],[42,44,46,48],[17,21,67,55],[15,27,65,53],[16,28,66,54]]
        l:L L I:=
         [[51,53,55,57],[52,54,56,58],[11,41,65,35],[17,47,63,33],[18,48,64,34]]
        b:L L I:=
         [[61,63,65,67],[62,64,66,68],[45,25,35,55],[47,27,37,57],[46,26,36,56]]
        llli2gp [f,r,u,d,l,b]

      mathieu11(l:L I):PERMGRP I ==
      -- permutations derived from the ATLAS
        l:=removeDuplicates l
        #l ~= 11 => error "Exactly 11 integers for mathieu11 needed !"
        a:L L I:=[[l.1,l.10],[l.2,l.8],[l.3,l.11],[l.5,l.7]]
        llli2gp [a,[[l.1,l.4,l.7,l.6],[l.2,l.11,l.10,l.9]]]

      mathieu11():PERMGRP I == mathieu11 li1n 11

      mathieu12(l:L I):PERMGRP I ==
      -- permutations derived from the ATLAS
        l:=removeDuplicates l
        #l ~= 12 => error "Exactly 12 integers for mathieu12 needed !"
        a:L L I:=
          [[l.1,l.2,l.3,l.4,l.5,l.6,l.7,l.8,l.9,l.10,l.11]]
        llli2gp [a,[[l.1,l.6,l.5,l.8,l.3,l.7,l.4,l.2,l.9,l.10],[l.11,l.12]]]

      mathieu12():PERMGRP I == mathieu12 li1n 12

      mathieu22(l:L I):PERMGRP I ==
      -- permutations derived from the ATLAS
        l:=removeDuplicates l
        #l ~= 22 => error "Exactly 22 integers for mathieu22 needed !"
        a:L L I:=[[l.1,l.2,l.4,l.8,l.16,l.9,l.18,l.13,l.3,l.6,l.12],   _
          [l.5,l.10,l.20,l.17,l.11,l.22,l.21,l.19,l.15,l.7,l.14]]
        b:L L I:= [[l.1,l.2,l.6,l.18],[l.3,l.15],[l.5,l.8,l.21,l.13],   _
          [l.7,l.9,l.20,l.12],[l.10,l.16],[l.11,l.19,l.14,l.22]]
        llli2gp [a,b]

      mathieu22():PERMGRP I == mathieu22 li1n 22

      mathieu23(l:L I):PERMGRP I ==
      -- permutations derived from the ATLAS
        l:=removeDuplicates l
        #l ~= 23 => error "Exactly 23 integers for mathieu23 needed !"
        a:L L I:= [[l.1,l.2,l.3,l.4,l.5,l.6,l.7,l.8,l.9,l.10,l.11,l.12,l.13,l.14,_
                   l.15,l.16,l.17,l.18,l.19,l.20,l.21,l.22,l.23]]
        b:L L I:= [[l.2,l.16,l.9,l.6,l.8],[l.3,l.12,l.13,l.18,l.4],              _
                   [l.7,l.17,l.10,l.11,l.22],[l.14,l.19,l.21,l.20,l.15]]
        llli2gp [a,b]

      mathieu23():PERMGRP I == mathieu23 li1n 23

      mathieu24(l:L I):PERMGRP I ==
      -- permutations derived from the ATLAS
        l:=removeDuplicates l
        #l ~= 24 => error "Exactly 24 integers for mathieu24 needed !"
        a:L L I:= [[l.1,l.16,l.10,l.22,l.24],[l.2,l.12,l.18,l.21,l.7],          _
                   [l.4,l.5,l.8,l.6,l.17],[l.9,l.11,l.13,l.19,l.15]]
        b:L L I:= [[l.1,l.22,l.13,l.14,l.6,l.20,l.3,l.21,l.8,l.11],[l.2,l.10],  _
                   [l.4,l.15,l.18,l.17,l.16,l.5,l.9,l.19,l.12,l.7],[l.23,l.24]]
        llli2gp [a,b]

      mathieu24():PERMGRP I == mathieu24 li1n 24

      janko2(l:L I):PERMGRP I ==
      -- permutations derived from the ATLAS
        l:=removeDuplicates l
        #l ~= 100 => error "Exactly 100 integers for janko2 needed !"
        a:L L I:=[                                                            _
                 [l.2,l.3,l.4,l.5,l.6,l.7,l.8],                               _
                 [l.9,l.10,l.11,l.12,l.13,l.14,l.15],                         _
                 [l.16,l.17,l.18,l.19,l.20,l.21,l.22],                        _
                 [l.23,l.24,l.25,l.26,l.27,l.28,l.29],                        _
                 [l.30,l.31,l.32,l.33,l.34,l.35,l.36],                        _
                 [l.37,l.38,l.39,l.40,l.41,l.42,l.43],                        _
                 [l.44,l.45,l.46,l.47,l.48,l.49,l.50],                        _
                 [l.51,l.52,l.53,l.54,l.55,l.56,l.57],                        _
                 [l.58,l.59,l.60,l.61,l.62,l.63,l.64],                        _
                 [l.65,l.66,l.67,l.68,l.69,l.70,l.71],                        _
                 [l.72,l.73,l.74,l.75,l.76,l.77,l.78],                        _
                 [l.79,l.80,l.81,l.82,l.83,l.84,l.85],                        _
                 [l.86,l.87,l.88,l.89,l.90,l.91,l.92],                        _
                 [l.93,l.94,l.95,l.96,l.97,l.98,l.99] ]
        b:L L I:=[
                [l.1,l.74,l.83,l.21,l.36,l.77,l.44,l.80,l.64,l.2,l.34,l.75,l.48,l.17,l.100],_
                [l.3,l.15,l.31,l.52,l.19,l.11,l.73,l.79,l.26,l.56,l.41,l.99,l.39,l.84,l.90],_
                [l.4,l.57,l.86,l.63,l.85,l.95,l.82,l.97,l.98,l.81,l.8,l.69,l.38,l.43,l.58],_
                [l.5,l.66,l.49,l.59,l.61],_
                [l.6,l.68,l.89,l.94,l.92,l.20,l.13,l.54,l.24,l.51,l.87,l.27,l.76,l.23,l.67],_
                [l.7,l.72,l.22,l.35,l.30,l.70,l.47,l.62,l.45,l.46,l.40,l.28,l.65,l.93,l.42],_
                [l.9,l.71,l.37,l.91,l.18,l.55,l.96,l.60,l.16,l.53,l.50,l.25,l.32,l.14,l.33],_
                [l.10,l.78,l.88,l.29,l.12] ]
        llli2gp [a,b]

      janko2():PERMGRP I == janko2 li1n 100

      abelianGroup(l:L PI):PERMGRP I ==
        gens:= nil()$(L L L I)
        element:I:= 1
        for n in l | n > 1 repeat
          gens:=cons( list [i for i in element..(element+n-1) ], gens )
          element:=element+n
        llli2gp
          #gens = 0 => [[[1]]]
          gens

      alternatingGroup(l:L I):PERMGRP I ==
        l:=removeDuplicates l
        #l = 0 =>
          error "Cannot construct alternating group on empty set"
        #l < 3 => llli2gp [[[l.1]]]
        #l = 3 => llli2gp [[[l.1,l.2,l.3]]]
        tmp:= [l.i for i in 3..(#l)]
        gens:L L L I:=[[tmp],[[l.1,l.2,l.3]]]
        odd?(#l) => llli2gp gens
        gens.1 := cons([l.1,l.2],gens.1)
        llli2gp gens

      alternatingGroup(n:PI):PERMGRP I == alternatingGroup li1n n

      symmetricGroup(l:L I):PERMGRP I ==
        l:=removeDuplicates l
        #l = 0 => error "Cannot construct symmetric group on empty set !"
        #l < 3 => llli2gp [[l]]
        llli2gp [[l],[[l.1,l.2]]]

      symmetricGroup(n:PI):PERMGRP I == symmetricGroup li1n n

      cyclicGroup(l:L I):PERMGRP I ==
        l:=removeDuplicates l
        #l = 0 => error "Cannot construct cyclic group on empty set"
        llli2gp [[l]]

      cyclicGroup(n:PI):PERMGRP I == cyclicGroup li1n n

      dihedralGroup(l:L I):PERMGRP I ==
        l:=removeDuplicates l
        #l < 3 => error "in dihedralGroup: Minimum of 3 elements needed !"
        tmp := [[l.i, l.(#l-i+1) ] for i in 1..(#l quo 2)]
        llli2gp [ [ l ], tmp ]

      dihedralGroup(n:PI):PERMGRP I ==
        n = 1 => symmetricGroup (2::PI)
        n = 2 => llli2gp [[[1,2]],[[3,4]]]
        dihedralGroup li1n n

