/usr/share/axiom-20170501/src/algebra/PERMGRP.spad

)abbrev domain PERMGRP PermutationGroup
++ Authors: G. Schneider, H. Gollan, J. Grabmeier
++ Date Created: 13 February 1987
++ Date Last Updated: 24 May 1991
++ References:
++ Sims71 Determining the conjugacy classes of a permutation group
++ Description: 
++ PermutationGroup implements permutation groups acting
++ on a set S, 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 Axiom 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) : SIG == CODE where
  S : SetCategory

  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 )
    -- REC holds orbit and Schreier vector
  REC2 ==> Record(order:NNI,sgset:L V NNI,_
             gpbase:L NNI,orbs:L REC,mp:L S,wd:L L NNI)
    -- REC2 holds extra information about group in representation
    --      to improve efficiency of some functions.
    --      See Rep below for more details.
  REC3 ==> Record(elt:V NNI,lst:L NNI)
    -- REC3 holds an element and a word
  REC4 ==> Record(bool:B,lst:L NNI)
    -- REC4 used by 'memberInternal' function to return internal
    -- membership testing

  SIG ==> SetCategory with

    coerce : % -> L PERM S
      ++ coerce(gp) returns the generators of the group gp.
      ++
      ++X x : PERM INT := [[1,3,5],[7,11,9]]

    coerce : L PERM S -> %
      ++ coerce(ls) coerces a list of permutations ls to the group
      ++ generated by this list.
      ++
      ++X y : PERM INT := [[3,5,7,9]]
      ++X z : PERM INT := [1,3,11]
      ++X g : PERMGRP INT := [ y , z ]

    generators : % -> L PERM S
      ++ generators(gp) returns the generators of the group gp.
      ++
      ++X S:List(Integer) := [1,2,3,4]
      ++X G := symmetricGroup(S)
      ++X generators(G)

    elt : (%,NNI) -> PERM S
      ++ elt(gp,i) returns the i-th generator of the group gp.
      ++
      ++X S:List(Integer) := [1,2,3,4]
      ++X G := symmetricGroup(S)
      ++X G.2

    random : (%,I) -> PERM S
      ++ random(gp,i) returns a random product of maximal i generators
      ++ of the group gp.
      ++
      ++X S:List(Integer) := [1,2,3,4]
      ++X G := symmetricGroup(S)
      ++X random(G,2)

    random : % -> PERM S
      ++ random(gp) returns a random product of maximal 20 generators
      ++ of the group gp.
      ++ Note: random(gp)=random(gp,20).
      ++
      ++X S:List(Integer) := [1,2,3,4]
      ++X G := symmetricGroup(S)
      ++X random(G)

    order : % -> NNI
      ++ order(gp) returns the order of the group gp.
      ++
      ++X x : PERM INT := [[1,3,5],[7,11,9]]
      ++X y : PERM INT := [[3,5,7,9]]
      ++X g : PERMGRP INT := [ x , y ]
      ++X order g

    degree : % -> NNI
      ++ degree(gp) returns the number of points moved by all permutations
      ++ of the group gp.
      ++
      ++X y : PERM INT := [[3,5,7,9]]
      ++X z : PERM INT := [1,3,11]
      ++X g : PERMGRP INT := [ y , z ]
      ++X degree g

    base : % -> L S
      ++ base(gp) returns a base for the group gp.
      ++
      ++X S:List(Integer) := [1,2,3,4]
      ++X G := symmetricGroup(S)
      ++X base(G)

    strongGenerators : % -> L PERM S
      ++ strongGenerators(gp) returns strong generators for
      ++ the group gp.
      ++
      ++X S:List(Integer) := [1,2,3,4]
      ++X G := symmetricGroup(S)
      ++X strongGenerators(G)

    wordsForStrongGenerators : % -> L L NNI
      ++ wordsForStrongGenerators(gp) returns the words for the strong
      ++ generators of the group gp in the original generators of
      ++ gp, represented by their indices in the list, given by
      ++ generators.

    permutationGroup : L PERM S -> %
      ++ permutationGroup(ls) coerces a list of permutations ls to
      ++ the group generated by this list.

    orbit : (%,S) -> FSET S
      ++ orbit(gp,el) returns the orbit of the element el under the
      ++ group gp, the set of all points gained by applying
      ++ each group element to el.
      ++
      ++X S:List(Integer) := [1,2,3,4]
      ++X G := symmetricGroup(S)
      ++X orbit(G,[1,2,3])
      ++X x : PERM INT := [[1,3,5],[7,11,9]]
      ++X y : PERM INT := [[3,5,7,9]]
      ++X g : PERMGRP INT := [ x , y ]
      ++X orbit(g, 3)

    orbit : (%,FSET S)-> FSET FSET S
      ++ orbit(gp,els) returns the orbit of the unordered
      ++ set els under the group gp.

    orbit : (%,L S) -> FSET L S
      ++ orbit(gp,ls) returns the orbit of the ordered
      ++ list ls under the group gp.
      ++ Note: return type is L L S temporarily because FSET L S has an error.
      ++
      ++X S:List(Integer) := [1,2,3,4]
      ++X G := symmetricGroup(S)
      ++X orbit(G,[1,2,3])

    orbits : % -> FSET FSET S
      ++ orbits(gp) returns the orbits of the group gp,
      ++ it partitions the (finite) of all moved points.
      ++
      ++X y : PERM INT := [[3,5,7,9]]
      ++X z : PERM INT := [1,3,11]
      ++X g : PERMGRP INT := [ y , z ]
      ++X orbits g

    member? : (PERM S, %)-> B
      ++ member?(pp,gp) answers the question, whether the
      ++ permutation pp is in the group gp or not.
      ++
      ++X x : PERM INT := [[1,3,5],[7,11,9]]
      ++X y : PERM INT := [[3,5,7,9]]
      ++X z : PERM INT := [1,3,11]
      ++X g : PERMGRP INT := [ x , z ]
      ++X member? ( y , g )

    wordInStrongGenerators : (PERM S, %)-> L NNI
      ++ wordInStrongGenerators(p,gp) returns the word for the
      ++ permutation p in the strong generators of the group gp,
      ++ represented by the indices of the list, given by strongGenerators.

    wordInGenerators : (PERM S, %)-> L NNI
      ++ wordInGenerators(p,gp) returns the word for the permutation p
      ++ in the original generators of the group gp,
      ++ represented by the indices of the list, given by generators.

    movedPoints : % -> FSET S
      ++ movedPoints(gp) returns the points moved by the group gp.
      ++
      ++X x : PERM INT := [[1,3,5],[7,11,9]]
      ++X z : PERM INT := [1,3,11]
      ++X g : PERMGRP INT := [ x , z ]
      ++X movedPoints g

    "<" : (%,%) -> B
      ++ gp1 < gp2 returns true if and only if gp1
      ++ is a proper subgroup of gp2.

    "<=" : (%,%) -> B
      ++ gp1 <= gp2 returns true if and only if gp1
      ++ is a subgroup of gp2.
      ++ Note: because of a bug in the parser you have to call this
      ++ function explicitly by gp1 <=$(PERMGRP S) gp2.
      -- (GILT DAS NOCH?)

    initializeGroupForWordProblem : % -> Void
      ++ initializeGroupForWordProblem(gp) initializes the group gp
      ++ for the word problem.
      ++ Notes: it calls the other function of this name with parameters
      ++ 0 and 1: 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
      ++ 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.

  CODE ==> add

    -- Representation of the instance:
    -- The 'gens' component completely defines the group as a list
    -- of permutations. This is set when the group is constructed.
    -- The information component allows some functions to be run
    -- more efficiently this data is created, when needed from gens.
    -- The parts of the information data are defined as follows:
    -- order  - Number of elements. Zero means that 'information'
    --          data has not yet been computed.
    -- sgset  - Strong Generators
    -- gpbase - sequence of points stabilised by the group.
    -- orbs   - Describes orbits of base point. The orb part is
    --          just list of point on the orbit. The Schreier vector
    --          (svc) part allows you to compute element of the group
    --          moving given point to base point of the orbit.
    --          -2 means point not in orbit,
    --          -1 means base point,
    --          positive value is index of strong generator moving
    --          given point closer to base point.
    --          This list of orbits tends to be in a certain order,
    --          (corresponding to the order of gpbase)
    --          that is, stabiliser of point 1 (if it exists) is first
    --          then the other stabilisers, then
    --          the final orbit may not stabilise any points.
    --          I don't know if this order is required or assumed
    --          by any functions.
    -- mp     - List of elements of S moved by some permutation
    --          (needed for mapping between permutations on S and
    --          internal representation)
    -- wd     - Gives representation of strong generators in terms
    --          of original generators

    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 ^entryLessZero repeat
        entry := schreierVector.(actelt.point)
        entryLessZero  := (entry < 0)
        if  ^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 ( eval ( gen , workList.j ) , newList )
          if ^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 maxLoops > 0 then
        numberOfLoops : I  := 1 + (random()$Integer rem maxLoops)
      else
        numberOfLoops : I := maxLoops
      while numberOfLoops > 0 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 ==
      support : FSET S :=  brace()   -- empty set !!
      for perm in group.gens repeat
        support := union(support, movedPoints perm)
      parts support

    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 ^ 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 degree > 0 repeat
        p := osvc.ppt
        p < 0 => return [ xelt , word ]
        x    := group.p
        xelt := times ( x , xelt )
        if wordProblem then word := append ( wordlist.p , word )
        ppt  := x.ppt

    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
      for i in number1..degree repeat
        ort := orbitWithSvc ( group , i )
        k   := ort.orb
        k1  := # k
        if k1 ^= 1 then leave
      gpsgs := nil()$(L V NNI)
      words2 := nil()$(L L NNI)
      gplength : NNI := #group
      for jj 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 j > 0 repeat
        -- find generators for the stabilizer
        ran := ranelt ( group , words , maxLoops )
        str := strip ( ran.elt , ort , group , words )
        el2 := str.elt
        if ^ testIdentity el2 then
          if ^ 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 ( maxLoops < 0 ) 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   := 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 )
      void()


    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 := eval ( 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
      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
        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 ( degree > 0 ) and noresult repeat
                  entry := s.(z.p)
                  if entry < 0 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 degree gen > 0 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
      void

    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 movedPoints 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   := eval ( 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

    movedPoints ( 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 numberOfLoops > 0 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    := movedPoints gp
      orbits := nil()$(L FSET S)
      while cardinality spp > 0 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 ( movedPoints gp1 , movedPoints gp2 ) => false
      not subgroup ( gp1 , gp2 ) => false
      order gp1 = order gp2 => false
      true

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

    gp1 = gp2 ==
      movedPoints gp1 ^= movedPoints 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 ]
      void

    initializeGroupForWordProblem ( gp ) ==
       initializeGroupForWordProblem ( gp , 0 , 1 )
axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / PERMGRP.spad
