gpresent.spad

)if false
\documentclass{article}
\usepackage{axiom}
\usepackage{url}
\begin{document}
\title{Group Presentations}
\author{Martin J Baker}
\maketitle
\begin{abstract}
This file implements group structures related to algebraic topology,
specifically its homotopy and homology.

There are two such structures in this file:
\begin{itemize}
\item GroupPresentation - Defines a group by its generators and
relations. Used to hold fundamental group (homotopy).
\item Homology - Intended to hold homology which is calculated using
integer Smith normal form. This is an abelian group.
\end{itemize}
I have put a fuller explanation of this code here:
\url{http://www.euclideanspace.com/prog/scratchpad/mycode/discrete/finiteGroup/presentation/}
\end{abstract}

\section{Introduction}
Group represented by its generators and relations.
Here we use it to hold homotopy group such as fundamental group.
for more documentation see:
\url{http://www.euclideanspace.com/prog/scratchpad/mycode/topology/}

Representation holds the group as a set of generators and a set of
relations
Each generator is a NNI
Each relation is a list of indexes to generators. Negative values
indicate the inverse of the generator. So if '1' represents the
generator 'A' then '-1' represents its inverse 'A^-1'.

Note that the use of negative indices to represent the inverse does
not imply an abelian group. This is just a convenient way to code
the representation and, in general, the group is not necessarily
abelian.
The notation uses the following conventions.
\begin{itemize}
\item If index is positive it represents a generator shown as an
alphabetic digit followed by number if required.
\item If index is zero it represents identity shown as 'e'.
\item If index is negative it represents an inverse generator.
The notation uses '-' to indicate an inverse generator (I know this is
an abuse of notation, because the group is not necessarily Abelian).
\end{itemize}
This convention is not ambiguous and I don't like the alternatives,
'-1' exponent uses too much space and upper case alpha for inverses
may not be clear to people who don't read documentation.

\section{Homotopy Group}
The homotopy group is finitely presented by generators and relations.
This representation of a group is not, in general, algorithmically
computable into other representations of a group.

We can therefore compute 'a' (not 'the') homotopy group for a given
simplicial complex. We may also be able to apply some simplifications
to this group. However, in the general case, we cannot determine if
this is the simplest representation or determine if two such groups
are isomophic (their corresponding simplicial complexes are
homeomorphic).

Despite these fundamental limits on what is theoretically possible
I still believe it is worthwhile to have the capability to generate
'a' homotopy group for a given structure.

\section{Conversion between Finitely Presented and Permutation Groups}

In FriCAS the main group functionality is in PermutationGroup so it
is useful to be able to convert to and from it. For conversions to
PermutationGroup we use a well known algorithm called the
Todd-Coxeter algorithm.

If you wish to see how the algorithm works (or does not work) try
calling the function, with trace set to 'true', like this:

\begin{verbatim}
(7) -> d3 := dihedralGroup(3)$GroupPresentation

   (7)  <a b |  a*a*a, b*b, a*b*a*b>
                            Type: GroupPresentation
(8) -> toPermutationIfCan(d3,true)

   addPoint: cannot deduce more so adding a point
   adding:1 at row:1
                  +2  2  2+ +0  0+ +2  0  2  0+
   relatorTables=[|       |,|    |,|          |]
                  +0  0  0+ +0  0+ +0  0  0  0+
   addPoint: cannot deduce more so adding a point
   adding:2 at row:1
                  +2  2  2+ +3  3+ +2  3  2  3+
                  |       | |    | |          |
   relatorTables=[|0  0  0|,|0  0|,|0  0  0  0|]
                  |       | |    | |          |
                  +0  0  0+ +0  0+ +0  0  0  0+
   inferFromRelations found a gap of 1 so deduction made
   gap is in relator table:[2,2]
   distance from start:1 from end:0
   row at start:1
   row to change:3 new value:1
   generator index=2 invm=false
   forwardSequence:[1,3] backwardSequence:[1,0]
   inferFromRelations genIn=2 gb=1
                  +2  2  2+ +3  3+ +2  3  2  3+
                  |       | |    | |          |
   relatorTables=[|0  0  0|,|0  0|,|0  0  0  0|]
                  |       | |    | |          |
                  +0  0  0+ +1  1+ +0  1  0  1+
   addPoint: cannot deduce more so adding a point
   adding:1 at row:2
                  +2  2  2+ +3  3+ +2  3  2  3+
                  |       | |    | |          |
                  |4  4  4| |0  0| |4  0  4  0|
   relatorTables=[|       |,|    |,|          |]
                  |0  0  0| |1  1| |0  1  0  1|
                  |       | |    | |          |
                  +0  0  0+ +0  0+ +0  0  0  0+
   inferFromRelations found a gap of 1 so deduction made
   gap is in relator table:[1,1,1]
   distance from start:2 from end:0
   row at start:1
   row to change:4 new value:1
   generator index=1 invm=false
   forwardSequence:[1,2,4] backwardSequence:[1,0,0]
   inferFromRelations genIn=1 gb=2
                  +2  2  2+ +3  3+ +2  3  2  3+
                  |       | |    | |          |
                  |4  4  4| |0  0| |4  0  4  0|
   relatorTables=[|       |,|    |,|          |]
                  |0  0  0| |1  1| |0  1  0  1|
                  |       | |    | |          |
                  +1  1  1+ +0  0+ +1  0  1  0+
   addPoint: cannot deduce more so adding a point
   adding:2 at row:2
                  +2  2  2+ +3  3+ +2  3  2  3+
                  |       | |    | |          |
                  |4  4  4| |5  5| |4  5  4  5|
                  |       | |    | |          |
   relatorTables=[|0  0  0|,|1  1|,|0  1  0  1|]
                  |       | |    | |          |
                  |1  1  1| |0  0| |1  0  1  0|
                  |       | |    | |          |
                  +0  0  0+ +0  0+ +0  0  0  0+
   inferFromRelations found a gap of 1 so deduction made
   gap is in relator table:[2,2]
   distance from start:1 from end:0
   row at start:2
   row to change:5 new value:2
   generator index=2 invm=false
   forwardSequence:[2,5] backwardSequence:[2,0]
   inferFromRelations genIn=2 gb=1
                  +2  2  2+ +3  3+ +2  3  2  3+
                  |       | |    | |          |
                  |4  4  4| |5  5| |4  5  4  5|
                  |       | |    | |          |
   relatorTables=[|0  0  0|,|1  1|,|0  1  0  1|]
                  |       | |    | |          |
                  |1  1  1| |0  0| |1  0  1  0|
                  |       | |    | |          |
                  +0  0  0+ +2  2+ +0  2  0  2+
   inferFromRelations found a gap of 1 so deduction made
   gap is in relator table:[1,2,1,2]
   distance from start:2 from end:1
   row at start:1
   row to change:5 new value:3
   generator index=1 invm=false
   forwardSequence:[1,2,5,0] backwardSequence:[1,3,0,0]
   inferFromRelations genIn=1 gb=2
                  +2  2  2+ +3  3+ +2  3  2  3+
                  |       | |    | |          |
                  |4  4  4| |5  5| |4  5  4  5|
                  |       | |    | |          |
   relatorTables=[|0  0  0|,|1  1|,|0  1  0  1|]
                  |       | |    | |          |
                  |1  1  1| |0  0| |1  0  1  0|
                  |       | |    | |          |
                  +3  3  3+ +2  2+ +3  2  3  2+
   addPoint: cannot deduce more so adding a point
   adding:1 at row:3
                  +2  2  2+ +3  3+ +2  3  2  3+
                  |       | |    | |          |
                  |4  4  4| |5  5| |4  5  4  5|
                  |       | |    | |          |
                  |6  6  6| |1  1| |6  1  6  1|
   relatorTables=[|       |,|    |,|          |]
                  |1  1  1| |0  0| |1  0  1  0|
                  |       | |    | |          |
                  |3  3  3| |2  2| |3  2  3  2|
                  |       | |    | |          |
                  +0  0  0+ +0  0+ +0  0  0  0+
   inferFromRelations found a gap of 1 so deduction made
   gap is in relator table:[1,1,1]
   distance from start:1 from end:1
   row at start:3
   row to change:6 new value:5
   generator index=1 invm=false
   forwardSequence:[3,6,0] backwardSequence:[3,5,0]
   inferFromRelations genIn=1 gb=1
                  +2  2  2+ +3  3+ +2  3  2  3+
                  |       | |    | |          |
                  |4  4  4| |5  5| |4  5  4  5|
                  |       | |    | |          |
                  |6  6  6| |1  1| |6  1  6  1|
   relatorTables=[|       |,|    |,|          |]
                  |1  1  1| |0  0| |1  0  1  0|
                  |       | |    | |          |
                  |3  3  3| |2  2| |3  2  3  2|
                  |       | |    | |          |
                  +5  5  5+ +0  0+ +5  0  5  0+
   inferFromRelations found a gap of 1 so deduction made
   gap is in relator table:[1,2,1,2]
   distance from start:1 from end:2
   row at start:2
   row to change:4 new value:6
   generator index=2 invm=false
   forwardSequence:[2,4,0,0] backwardSequence:[2,5,6,0]
   inferFromRelations genIn=2 gb=1
                  +2  2  2+ +3  3+ +2  3  2  3+
                  |       | |    | |          |
                  |4  4  4| |5  5| |4  5  4  5|
                  |       | |    | |          |
                  |6  6  6| |1  1| |6  1  6  1|
   relatorTables=[|       |,|    |,|          |]
                  |1  1  1| |6  6| |1  6  1  6|
                  |       | |    | |          |
                  |3  3  3| |2  2| |3  2  3  2|
                  |       | |    | |          |
                  +5  5  5+ +0  0+ +5  0  5  0+
   inferFromRelations found a gap of 1 so deduction made
   gap is in relator table:[2,2]
   distance from start:1 from end:0
   row at start:4
   row to change:6 new value:4
   generator index=2 invm=false
   forwardSequence:[4,6] backwardSequence:[4,0]
   inferFromRelations genIn=2 gb=1
                  +2  2  2+ +3  3+ +2  3  2  3+
                  |       | |    | |          |
                  |4  4  4| |5  5| |4  5  4  5|
                  |       | |    | |          |
                  |6  6  6| |1  1| |6  1  6  1|
   relatorTables=[|       |,|    |,|          |]
                  |1  1  1| |6  6| |1  6  1  6|
                  |       | |    | |          |
                  |3  3  3| |2  2| |3  2  3  2|
                  |       | |    | |          |
                  +5  5  5+ +4  4+ +5  4  5  4+

   (8)  <(1 2 4)(3 6 5),(1 3)(2 5)(4 6)>
                    Type: Union(PermutationGroup(Integer),...)
\end{verbatim}

This will display the relatorTables at each stage and the deductions
being made from the tables.

The conversions in both directions can be improved by implementing
better simplifications but since there is no canonical form the
simplifications can never be perfect. The Todd-Coxeter does not yet
detect and remove 'coincidences', that is, duplicated points. So
this in the next improvement that needs to be made.

\section{Direct Product of Groups}

If G = < Sg | Rg > and H = < Sh | Rh >
Then G*H = < Sg,Sh | Rg,Rh,Rp >

where:
G*H is the direct product of G and H.
Rp is a set of relations specifying that each element of Sg
anti-commutes with each element of Sh.
See:
\url{https://en.wikipedia.org/wiki/Direct_product_of_groups}

This assumes that Sg and Sh are disjoint. In order to assure this the
generators will be renumbered before doing the product.

\section{Quotients}

We can generate a quotient by removing a generator (more precisely
adding relation equating this generator to neutral element) or
more generaly by adding a relation.

A quotient is an example of a surjective mapping
(epimorphism) on the group. It is easy to work with functions
between finitely presented groups, all we have to do is map
the generators and check that the relations still equal 1.

\section{Simplify Finitely Presented Groups}

There may not be a simplest form but it is possible to do some
simplifications.

In order to try to simplify a finitely generated presentation
to produce simpler but isomorphic groups we can apply certain
transformations or automorphisms (isomorphisms of the group
back to itself).

For example:
\begin{itemize}
\item Remove all zero terms in relations.
\item If a relation consists of a single generator then remove
  that generator.
\item If a relation consists of a pair of generators then make
  the second generator the inverse of the first.
\item If a generator is adjacent to its inverse then cancel
  them out.
\item Remove duplicate relations.
\item Substitute one relation in another.
\end{itemize}

These automorphisms were studied and categorised by Tietze
and Nielsen.

\subsection{Tietze Transformations}

\begin{table}[]
\label{Tietze transformations are of 4 kinds}
\begin{tabular}{lll}
 \ Kind \ Examples \\
T1 \ Add a relation \ < A | A^3 > -> < A | A^3 , A^6 > \\
T2 \ Remove a relation \ for example we can reverse the above
< A | A^3 , A^6> -> < A | A^3> \\
T3 \ Add a generator \ < A | A^3 > -> <A , B | A^3, B = A^2 > \\
T4 \ Remove a generator \ for example we can reverse the above
<A , B | A^3, B = A^2 > -> < A | A^3 > \\
\end{tabular}
\end{table}
We are interested in simplifying so we are mostly interested
in T2 and T4.
\subsection{T2}

T2 allows us to remove a relation but not generators. This
happens when a relation is redundant, that is it contains no
additional information than is already contained in the
other relations.

This happens, for example, where:

One relation is a multiple of another - In this case we can
remove the highest multiple but not the lowest.

One relation is the inverse of another - In this case we can
remove any one, but not both.

We can also simplify relations, for example, where an element
and its inverse are next to each other they can be
canceled out and removed.
\subsection{T4}

T4 allows us to remove a generator and corresponding relations.
\subsection{Nielsen Transformations}

The following transformations on a finitely generated free group
produce isomorphic groups.

\begin{itemize}
\item Switch A and B
\item Cyclically permute A, B, ... to B, ..., A.
\item Replace A with A^(-1)
\item Replace A with A*B
\item Substitute one relation in another
\end{itemize}

Most of these rules are self explanatory except substitute which
does the following:

If a generator is contained in exactly 2 relations then we may be
able to substitute one relation in another and remove that generator.
If, in at least one of these relations, the generator is contained
only once then we can move it to one side of the equation and
substitute it in the other relation.

This is done by a local function TTSubstitute.

Here is an example of its use from SimplicialComplex without
substitution:
\begin{verbatim}
(1) -> tS := torusSurface()$SimplicialComplexFactory

     (1)
           (1,2,3)
           (2,3,5)
           (2,4,5)
           (2,4,7)
           (1,2,6)
           (2,6,7)
           (3,4,6)
           (3,5,6)
           (3,4,7)
           (1,3,7)
           (1,4,5)
           (1,4,6)
           (5,6,7)
           (1,5,7)
               Type: FiniteSimplicialComplex(VertexSetAbstract)

(2) -> fundamentalGroup(tS)

     (2)  <o t w |   o*w*t,  o*t*w>
                                        Type: GroupPresentation
\end{verbatim}
This needs to be further simplified,

<o t w |   -o = w*t,  -o = t*w>

Substituting for -o we have:

w*t = t*w

That is, two edges that commute.

Moving everything to one side of the equation we have:
t * w * -t * -w
\begin{verbatim}
(2) -> fundamentalGroup(tS)

   (2)  <s v |  -s*v*s*-v>
                                        Type: GroupPresentation
\end{verbatim}
This is the same as above with v=w and s= -t.

\section{Testing and Validating this Code}

Some functions are very difficult to test, for example in
the SimplicialComplex code here:
\url{http://www.euclideanspace.com/prog/scratchpad/mycode/topology/simplex/}
and here:
\url{http://www.euclideanspace.com/prog/scratchpad/mycode/topology/delta/}
there are functions such as fundamentalGroup which output a
GroupPresentation. The reason for the difficulty is that they
do not have a canonical form, that is, there may
be more than one correct result and none of them are better
than the others and there is no general algorithm for testing
if they are equal.

So some change to the code may change the result but the result
may be just as correct as the other result. So testing that
fundamentalGroup generates a given output for a given input
is not a useful test for correctness.

I think that all we can do in this situation is to test
fundamentalGroup with very simple inputs such as a
topological sphere. This should always produce an empty
presentation.
)endif

)abbrev domain GROUPP GroupPresentation
++ Author: Martin Baker
++ Description:
++   Group represented by its generators and relations.
++   Here we use it to hold homotopy group such as fundamental group.
++   for more documentation see:
++   http://www.euclideanspace.com/prog/scratchpad/mycode/discrete/finiteGroup/presentation/
++ Date Created: Jan 2016
++ Basic Operations:
++ Related packages:
++ Related categories:
++ Related Domains: PermutationGroup
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:

GroupPresentation() : Exports == Impl where
  NNI ==> NonNegativeInteger
  PI ==> PositiveInteger
  x<<y ==> hconcat(x::OutputForm, y::OutputForm)
  GENMAP ==> List(Record(OldGen : NNI, NewGen : NNI))
  Exports ==> SetCategory() with
    groupPresentation : (v : List(NNI), rels1 : List(List(Integer))) -> %
      ++ construct from generators and relations
    groupPresentation : (v : List(NNI)) -> %
      ++ construct free group with generators but no relations
    groupPresentation : () -> %
      ++ construct trivial group with no generators or relations
    simplify : (s : %) -> %
      ++ simplify(s) tries to simplify s.
      ++ There may not be a simplest form but it is possible to
      ++ do some simplifications as follows:
      ++ 1. Remove all zero terms in relations.
      ++ 2. If a relation consists of a single generator then remove
      ++    that generator.
      ++ 3. If a relation consists of a pair of generators then make the
      ++    second generator the inverse of the first.
      ++ 4. If a generator is adjacent to its inverse then cancel them out.
      ++ 5. Remove duplicate relations.
      ++ 6. Substitute one relation in another.
    simplify : (s : %, trace : Boolean) -> %
      ++ simplify with option to trace
    refactor : (a : %) -> %
      ++ actual value of generators is not important, it is only important that
      ++ they correspond to the appropriate entries in the relations.
      ++ Therefore we can refactor the generators without changing the
      ++ group represented.
    quotient : (a : %, remgen : List(NNI)) -> %
      ++ take quotient by removing generators specified by remgen
    quotient : (a : %, addrel : List(List(Integer))) -> %
      ++ take quotient by adding relations specified by addrel
    directProduct : (a : %, b : %) -> %
      ++ directProduct of two groups
    cyclicGroup : (n : PI) -> %
      ++ cyclicGroup(n) constructs the cyclic group of order n acting
      ++ on the integers 1, ..., n.
    dihedralGroup : (n : PI) -> %
      ++ dihedralGroup(n) constructs the dihedral group of order 2n
      ++ acting on integers 1, ..., N.
    symmetricGroup : (n : PI) -> %
      ++ symmetricGroup(n) constructs the symmetric group of order n-1.
      ++ Note: generates all possible relations may not be minimal.
    toPermutationIfCan : (a : %) -> Union(PermutationGroup Integer, "failed")
      ++ convert to permutation group return "failed" for infinite groups.
      ++ For more information about the algorithm see:
      ++ \url{http://www.euclideanspace.com/prog/scratchpad/mycode/discrete/finiteGroup/pres2perm/}
    toPermutationIfCan : (a : %, trace : Boolean)
                          -> Union(PermutationGroup Integer, "failed")
      ++ convert to permutation group return "failed" for infinite groups.
      ++ For more information about the algorithm see:
      ++ \url{http://www.euclideanspace.com/prog/scratchpad/mycode/discrete/finiteGroup/pres2perm/}
    toPermutationIfCan : (a : %, sg : List(List(Integer)), trace : Boolean)
                          -> Union(PermutationGroup Integer, "failed")
      ++ toPermutationIfCan(a, sg, trace) returns permutation
      ++ representation of a on cosets of subgroup generate by sg
      ++ or "failed" if computation exceed resource limit.
      ++ trace activates debugging printouts.

  Impl ==> add

   -- Representation holds the group as a set of generators and a set of
   -- relations.
   -- Each generator is a NNI
   -- Each relation is a list of indexes to generators.
   -- if index is positive it represents a generator output as an
   -- alphabetic digit followed by number if required.
   -- if index is zero it represents identity output as 'e'.
   -- if index is negative it represents an inverse generator.

   Rep := Record(gens : PrimitiveArray(NNI), rels : List(List(Integer)))

   -- construct from generators and relations
   groupPresentation(gens1 : List(NNI), rels1 : List(List(Integer))) : % ==
       -- print("  groupPresentation construct (" << gens1 << ", "
       --       << rels1<< ")")
       g := construct(gens1)$PrimitiveArray(NNI)
       --print("groupPresentation constuct : " << rels1)
       -- remove empty relations since this simplifies equality function
       [g, [r for r in rels1 | not(empty?(r))]]

   -- construct free group with generators but no relations
   groupPresentation(gens1 : List(NNI)) : % ==
       -- print("  groupPresentation construct (" << gens1 << ", "
       --       << rels1<< ")")
       g := construct(gens1)$PrimitiveArray(NNI)
       rels2 := []$List(List(Integer))
       [g, rels2]

   -- construct trivial group with no generators or relations
   groupPresentation() : % ==
       gens1 := []$List(NNI)
       rels1 := []$List(List(Integer))
       groupPresentation(gens1, rels1)

   -- Local function used by refactor to map a given generator in a relation.
   mapGen(a : Integer, ms : GENMAP) : Integer ==
       for m  in ms repeat
           if abs(a) = m.OldGen then return m.NewGen
           if abs(a) = -m.OldGen then return -m.NewGen
       error concat(["cant map ", string(a), " in refactor"])
       a

   -- Actual value of generators is not important, it is only important that
   -- they correspond to the appropriate entries in the relations.
   -- Therefore we can refactor the generators (to be ascending integers
   -- starting as 1) without changing the group represented.
   refactor(a : %) : % ==
       -- first generate a map from existing generators to new generators
       gms : GENMAP := empty()$GENMAP
       for g in entries(a.gens) for gn in 1..(#(a.gens)) repeat
           gm : Record(OldGen : NNI, NewGen : NNI) := [g, gn]
           gms := concat(gms, gm)
       -- now use this map to change elements of relations
       rels1 := []$List(List(Integer))
       for rel in a.rels repeat
           newRel : List(Integer) := empty()$List(Integer)
           for ele in rel repeat
               newEle : Integer := mapGen(ele, gms)
               newRel := concat(newRel, newEle)
           rels1 := concat(rels1, newRel)
       gens1 : List(NNI) := [gn for gn in 1..(#(a.gens))]
       groupPresentation(gens1, rels1)

   -- Isomorphism is the most useful level of 'equality' for
   -- groups but unfortunately this is not computable in
   -- the general case for presentations.
   -- Although exact equality is less useful it is still useful to
   -- compare very simple presentations in the validation code which
   -- is useful to give some level of confidence that the
   -- correct presentation was generated.
   -- TODO result can be dependent on initial generator order, for
   -- example <a, b | a*a, b*b*b> = <b, a | a*a, b*b*b> would be false
   -- should really check all permutations of generators and return
   -- true if any of them gives equality.
   _=(a : %, b : %) : Boolean ==
       ar : % := refactor(a)
       br : % := refactor(b)
       ags : List(NNI) := entries(ar.gens)
       bgs : List(NNI) := entries(br.gens)
       if set(ags)$Set(NNI) ~= set(bgs)$Set(NNI) then return false
       ars : List(List(Integer)) := entries(ar.rels)
       brs : List(List(Integer)) := entries(br.rels)
       set(ars)$Set(List(Integer)) = set(brs)$Set(List(Integer))

   -- Display generators as alphabetic digits
   -- Local function used by coerce to OutputForm and other functions.
   -- if i2 is positive it represents a generator shown as an
   -- alphabetic digit followed by number if required.
   -- if i2 is zero it represents identity shown as 'e'.
   -- if i2 is negative it represents an inverse generator.
   -- The notation uses '-' to indicate an inverse generator (I know
   -- this is an abuse of notation, because the group is not
   -- necessarily Abelian). I don't like the alternatives,
   -- '-1' exponent uses too much space and upper case for inverses
   -- may not be clear to people who don't read documentation.
   outputGen(i2 : Integer) : OutputForm ==
       (suffix, i) := divide(abs(i2), 25)
       letters : String := "eabcdfghijklmnopqrstuvwxyz"
       n : OutputForm := (letters(i + 1))::OutputForm
       -- print("  groupPresentation outputGen(" << i2 <<
       --       ") gives " << n)
       if suffix > 0 then n := hconcat(n, outputForm(suffix + 1))
       if i2 < 0 then return hconcat(message"-", n)
       n

   -- display a relation using alphabetic digits
   outputRel(r : List(Integer)) : OutputForm ==
       eleout : OutputForm := message("")
       seperator : OutputForm := message(" ")
       for ele in r repeat
           newterm : OutputForm := outputGen(ele)
           eleout := hconcat([eleout, seperator, newterm])$OutputForm
           seperator := message("*")
       eleout

   -- display a list of relations using alphabetic digits
   outputRelList(i2 : List(List(Integer))) : OutputForm ==
       rels1 : List(OutputForm) := []$List(OutputForm)
       for r in i2 repeat
           rels1 := concat(rels1,outputRel(r))
       if #rels1 > 0 then return commaSeparate(rels1)
       message(" ")

   -- display a list of generators using alphabetic digits
   outputGenList(ps : List(NNI)) : OutputForm ==
       gens1 : List(OutputForm) := []$List(OutputForm)
       for p in ps repeat
           gens1 := concat(gens1, outputGen(p::Integer))
       if #gens1 > 0 then return blankSeparate(gens1)
       message(" ")

   -- local function to return indexes to each relation containing a given
   -- generator.
   indexesOfRelUsingGen(s : %, gen : NNI) : List(NNI) ==
       res : List(NNI) := []
       r : List(List(Integer)) := s.rels
       for rel in r for reln in 1..(#r) repeat
           if member?(gen::Integer,rel) then res := concat(res,reln)
           if gen > 0 and member?(-(gen::Integer),rel)
               then res := concat(res,reln)
       res

   -- local function to remove generator 'val' from generators
   removeGen(gens1 : PrimitiveArray(NNI), val : NNI) : PrimitiveArray(NNI) ==
       remove(val, gens1)

   -- local function to remove generator 'val' from relations
   removeGen2(rels1 : List(List(Integer)), val : NNI) : List(List(Integer)) ==
       [remove(-val, remove(val::Integer, rel)) for rel in rels1]

   -- local function to replace generator 'val1' with 'val2'
   -- in relations
   replaceGen(rels1 : List(List(Integer)), val1 : NNI, val2 : Integer
             ) : List(List(Integer)) ==
       --print("  groupPresentation replaceGen=" << rels1 << _
       --      " val1=" << val1 << " val2=" << val2)
       rels2 := []$List(List(Integer))
       for rel in rels1 repeat
           rel2 := []$List(Integer)
           for ele in rel repeat
               e : Integer := abs(ele)
               if e = val1 then e := val2
               if ele < 0 then e := -e
               rel2 := concat(rel2, e)
           rels2 := concat(rels2, rel2)
       rels2

   -- Tietze Transformation to remove a generator that is equal to
   -- the identity element. That is there is a relation containing only one
   -- generator.
   -- This procedure searches for a single element relation, if found, it
   -- removes the corresponding generator and also removes it from
   -- any relations containing it.
   -- This procedure only removes one generator, if there are several
   -- such relations then this procedure needs to be called several times.
   -- This is a local function used by simplify.
   TTRemoveGeneratorIfIdentity(s : %, trace : Boolean) : % ==
       gens1 : PrimitiveArray(NNI) := s.gens
       rels1 : List(List(Integer)) := s.rels
       toBeRemoved : NNI := 0
       for rel in rels1 repeat
           if #rel = 1 and toBeRemoved = 0 then
               toBeRemoved := qcoerce(abs(first(rel)))
       if toBeRemoved = 0 then return s
       if trace then
           print hconcat([message("simplify: generator '"), _
               outputGen(toBeRemoved), _
               message("' is identity so remove it")])
       gens1 := removeGen(gens1, toBeRemoved)
       rels1 := removeGen2(rels1, toBeRemoved)
       if trace then print outputRelList(rels1)
       [gens1, rels1]

   -- Tietze Transformation to rename a generator.
   -- If a relation consists of a pair of generators then make the
   -- second generator the inverse of the first.
   -- This procedure searches for a two element relation, if found, it
   -- replaces the second element with the inverse of the first.
   -- This procedure only replaces one generator, if there are several
   -- such relations then this procedure needs to be called several times.
   -- This is a local function used by simplify.
   TTRenameGenerator(s : %, trace : Boolean) : % ==
       gens1 : PrimitiveArray(NNI) := s.gens
       rels1 : List(List(Integer)) := s.rels
       replaceFrom : NNI := 0
       replaceTo : Integer := 0
       for rel in rels1 repeat
           if #rel = 2 and replaceFrom = 0 then
               replaceTo := second(rel)
               replaceFrom := qcoerce(abs(first(rel)))
               if first(rel) > 0 then replaceTo := -replaceTo
               -- don't replace an element with itself or its inverse
               if replaceFrom = abs(replaceTo) then replaceFrom := 0
       if replaceFrom=0 then return s
       if trace then
           print hconcat([message("simplify: generator '"), _
               outputGen(replaceFrom), _
               message("' is replaced by '"), _
               outputGen(replaceTo), _
               message("'")])
       gens1 := removeGen(gens1, replaceFrom)
       rels1 := replaceGen(rels1, replaceFrom, replaceTo)
       if trace then print outputRelList(rels1)
       [gens1, rels1]

   -- This is a local function used by simplify.
   TTRemoveEmpty(s : %, trace : Boolean) : % ==
       gens1 : PrimitiveArray(NNI) := s.gens
       rels1 : List(List(Integer)) := s.rels
       rels2 : List(List(Integer)) := empty()$List(List(Integer))
       for rel in rels1 repeat
           --print("  groupPresentation simplify rel=" << rel)
           if not empty?(rel) then
               rels2 := concat(rels2, rel)
       [gens1, rels2]

   -- This is a local function used by simplify.
   TTRemoveZero(s : %, trace : Boolean) : % ==
       gens1 : PrimitiveArray(NNI) := s.gens
       rels1 : List(List(Integer)) := s.rels
       gens1 := removeGen(gens1, 0)
       rels1 := removeGen2(rels1, 0)
       [gens1, rels1]


   -- This is a local function used by simplify.
   TTRemoveEleTimesInverse(s : %, trace : Boolean) : % ==
       gens1 : PrimitiveArray(NNI) := s.gens
       rels1 : List(List(Integer)) := s.rels
       --print("TTRemoveEleTimesInverse relations in =" << rels1)
       rels2 : List(List(Integer)) := empty()$List(List(Integer))
       changed : Boolean := false
       for rel in rels1 repeat
           --print("TTRemoveEleTimesInverse rel=" << rel)
           rel2 : List(Integer) := empty()$List(Integer)
           lastele : Integer := 0
           for ele in rel repeat
               if abs(ele) = abs(lastele) and sign(ele) ~= sign(lastele) then
                   if trace then print hconcat([_
                       message("simplify: generator '"), _
                       outputGen(ele), _
                       message("' is adjacent to its inverse")])
                   changed := true
                   lastele := 0
               else
                   if lastele ~= 0 then rel2 := concat(rel2, lastele)
                   lastele := ele
           if lastele ~= 0 then rel2 := concat(rel2, lastele)
           if not empty?(rel2) then rels2 := concat(rels2, rel2)
       if trace and changed then print outputRelList(rels2)
       [gens1, rels2]

   -- local function to invert relation. Used by TTSubstitute,
   -- TTMinimiseInverses and relationEquivalent.
   -- We invert each element and then reverse the order.
   -- A bit like De Morgan's laws
   invertRelation(relationIn : List(Integer)) : List(Integer) ==
       relationOut := []$List(Integer)
       for ele in relationIn repeat
           relationOut := concat(-ele, relationOut)
       relationOut

   -- local function to cycle relation. Used by relationEquivalent.
   -- The effect of a relation is not changed by cycling
   cycleRelation(relationIn : List(Integer)) : List(Integer) ==
       relationOut : List(Integer) := concat(relationIn.rest,relationIn.first)
       --print(message "cycleRelation " << relationIn << message " to " << relationOut)
       relationOut

   -- Local function to test equivalence of two relations.
   -- Used by TTRemoveDuplicateRelation.
   -- Relations are considered equivalent if they are identical or
   -- if they are the same after being cycled or if they are the
   -- same after being inverted.
   relationEquivalent(relA : List(Integer),relB : List(Integer)) : Boolean ==
       -- first filter out cases where relations are different lengths
       if #relA ~= #relB then return false
       -- test for equality
       if relA = relB then return true
       -- test for equality with inverted.
       if relA = invertRelation(relB) then return true
       -- test for equality with cycle
       relBCycle : List(Integer) := copy relB
       for n in 1..(#relA) repeat
           relBCycle := cycleRelation(relBCycle)
           if relA = relBCycle then return true
           if relA = invertRelation(relBCycle) then return true
       false

   -- This is a local function used by simplify.
   -- It looks for and removes any duplicated relations.
   -- Relations are considered duplicates if they are identical or
   -- if they are the same after being cycled or if they are the
   -- same after being inverted.
   TTRemoveDuplicateRelation(s : %, trace : Boolean) : % ==
       gens1 : PrimitiveArray(NNI) := s.gens
       rels1 : List(List(Integer)) := s.rels
       rels2 := []$List(List(Integer))
       --print(message "TTRemoveDuplicateRelation =" << rels1)
       for rela in rels1 for nrela in 1..(#rels1) repeat
           -- include relation
           include : Boolean := true
           for relb in rels1 for nrelb in 1..(#rels1) repeat
               if nrela > nrelb then
                   if relationEquivalent(rela,relb) then
                       include : Boolean := false
                       if trace then
                           m ==> "TTRemoveDuplicateRelation duplicate found "
                           print(message m << rela << message "=" << relb)
           if include then
               rels2 := concat(rels2, rela)
       [gens1, rels2]

   -- This is a local function used by simplify.
   -- If a relation contains more inverted elements that non-inverted
   -- elements then it is easier to read if we invert all the terms.
   TTMinimiseInverses(s : %, trace : Boolean) : % ==
       gens1 : PrimitiveArray(NNI) := s.gens
       rels1 : List(List(Integer)) := s.rels
       rels2 := []$List(List(Integer))
       for rel in rels1 repeat
           numInverts : NNI := 0
           numNonInverts : NNI := 0
           for ele in rel repeat
               if ele < 0 then
                   numInverts := numInverts + 1
               else
                   numNonInverts := numNonInverts + 1
           if numInverts > numNonInverts then
               rels2 := concat(rels2, invertRelation(rel))
           else
               rels2 := concat(rels2, rel)
       [gens1, rels2]

   -- This is a local function used by TTSubstitute.
   -- Counts the number of times a generator (or its inverse) occurs
   -- in a relation.
   generatorOccurrences(rel : List(Integer),gen : NNI) : NNI ==
       res : NNI := 0
       for g in rel repeat
           if g = gen then res := res + 1
           if gen > 0 and g = -gen then res := res + 1
       res

   -- local function to remove relations containing given generator
   removeRelations(rels1 : List(List(Integer)), val : NNI
                  ) : List(List(Integer)) ==
       res : List(List(Integer)) := []$List(List(Integer))
       for rel in rels1 repeat
           if (not member?(val,rel)) and (not member?(-val,rel))then
               res := concat(res,rel)
       res

   -- This is a local function used by simplify.
   -- If a generator is contained in exactly 2 relations then we may be
   -- able to substitute one relation in another and remove that generator.
   -- If, in at least one of these relations, the generator is contained
   -- only once then we can move it to one side of the equation and
   -- substitute it in the other relation.
   TTSubstitute(s : %, trace : Boolean) : % ==
       gs : List(NNI) := entries(s.gens)
       rs : List(List(Integer)) := s.rels
       r1 : List(Integer) := []
       r2 : List(Integer) := []
       n1 : NNI := 0
       n2 : NNI := 0
       genToBeRemoved : NNI := 0
       for g in gs repeat
           indexes : List(NNI) := indexesOfRelUsingGen(s, g)
           if #indexes = 2 and genToBeRemoved=0 then
               -- we have a candidate for substitution but, to be
               -- sure generator must occur once in a relation
               genToBeRemoved := g
               r1 := rs.(indexes.1)
               r2 := rs.(indexes.2)
               n1 := generatorOccurrences(r1,g)
               n2 := generatorOccurrences(r2,g)
               if n1 ~= 1::NNI then
                   -- swap first and second relations
                   r3 : List(Integer) := r1 ; n3 : NNI := n1
                   r1 := r2 ; n1 := n2
                   r2 := r3 ; n2 := n3
               if n1 ~= 1::NNI then
                   genToBeRemoved := 0
       if n1 ~= 1::NNI then return s
       -- If we have got to this point then we know a substitution
       -- is possible.
       if trace then
           print(message("simplify: TTSubstitute (") << s << message(")"))
           print(message("genToBeRemoved=") << outputGen(genToBeRemoved) << _
                 message(" r1=") << outputRel(r1) <<
                 message(" r2=") << outputRel(r2))
           print(message("n1=") << n1 << message(" n2=") << n2)
       restr : List(Integer) := r1
       prer  : List(Integer) := []
       found : Boolean := false
       genInverted : Boolean := false
       while (not empty?(restr)) and (not found)repeat
           x : Integer := first(restr)
           restr := rest(restr)
           if x=genToBeRemoved or x= -genToBeRemoved
               then
                   found := true
                   if x<0 then genInverted := true
               else prer := concat(prer,x)
       postr  : List(Integer) := []
       while not empty?(restr) repeat
           x : Integer := first(restr)
           restr := rest(restr)
           postr := concat(postr,x)
       replacement := concat(invertRelation(prer),invertRelation(postr))
       -- now substitute replacement for genToBeRemoved in r2
       if trace then
           print(message("we will substitute ") << outputRel(replacement) <<
                 message(" for ") << outputGen(genToBeRemoved) <<
                 message(" in ") << outputRel(r2))
       newRel : List(Integer) := []
       for x in r2 repeat
           if abs(x)=abs(genToBeRemoved)
             then
                 if genInverted
                     then newRel := concat(newRel,invertRelation(replacement))
                     else newRel := concat(newRel,replacement)
             else newRel := concat(newRel,x)
       if trace then print(message("newRel=") << outputRel(newRel))
       gens1 : PrimitiveArray(NNI) := s.gens
       rels1 : List(List(Integer)) := s.rels
       gens1 := removeGen(gens1, genToBeRemoved)
       rels1 := removeRelations(rels1, genToBeRemoved)
       rels1 := concat(rels1,newRel)
       if trace then print(message("gens=") << outputGenList(entries(gens1))
                           << message(" rels=") << outputRelList(rels1))
       [gens1, rels1]

   -- true if 'a' is simpler than 'b'.
   -- There may not be an absolute measure of whether one presentation
   -- is simpler than another but this procedure is used only in specific
   -- circumstances, that is where we have attempted to simplify the
   -- presentation and we want to test if it is actually simpler.
   -- We do this by testing if the number of generators or relations has
   -- reduced or if the complexity of the relations has reduced.
   -- This is a local function used by simplify.
   isSimpler?(a : %, b : %) : Boolean ==
       gensa : PrimitiveArray(NNI) := a.gens
       relsa : List(List(Integer)) := a.rels
       gensb : PrimitiveArray(NNI) := b.gens
       relsb : List(List(Integer)) := b.rels
       if #gensa < #gensb then return true
       if #relsa < #relsb then return true
       relationCompleityA : NNI := 0
       for rel in relsa repeat
           relationCompleityA := relationCompleityA + #rel
       relationCompleityB : NNI := 0
       for rel in relsb repeat
           relationCompleityB := relationCompleityB + #rel
       if relationCompleityA < relationCompleityB then return true
       false

   simplify(s : %) : % ==
       simplify(s, false)

   -- There may not be a simplest form but it is possible to
   -- do some simplifications as follows:
   -- 1) remove all zero terms in relations
   -- 2) if a relation consists of a single generator then remove that
   --    generator
   -- 3) if a relation consists of a pair of generators then make the
   --    second generator the inverse of the first
   -- 4) if a generator is adjacent to its inverse then cancel them out.
   -- 5) remove duplicate relations
   -- 6) Substitute one relation in another
   simplify(s : %, trace : Boolean) : % ==
       if trace then
           print(message("simplify(") << s << message(")"))
       res : % := s
       lastpass : % := s
       level : NNI := 0
       loopBreaker : NNI := 0
       while loopBreaker < 10000 repeat
           loopBreaker := loopBreaker + 1
           if level=0 then res := TTRemoveEmpty(res, trace)
           if level=1 then res := TTRemoveZero(res, trace)
           if level=2 then res := TTRemoveGeneratorIfIdentity(res, trace)
           if level=3 then res := TTRenameGenerator(res, trace)
           if level=4 then res := TTRemoveEleTimesInverse(res, trace)
           if level=5 then res := TTRemoveDuplicateRelation(res, trace)
           if level=6 then res := TTSubstitute(res, trace)
           if level=7 then return TTMinimiseInverses(res, trace)
           if isSimpler?(res, lastpass)
               then level := 0
               else level := level + 1
           if trace then
               print(message(" level=") << level << _
                     message(" loop=") << loopBreaker << _
                     message(" res=") << res)
           lastpass := res
       print(message("simplify excessive time - loop suspected") << s)
       res

   -- local function used by directProduct to offset indexes
   offsetIndexes(a : %, offset : NNI) : % ==
       ga : List(NNI) := entries(a.gens)
       ra : List(List(Integer)) := a.rels
       gb := [x+offset for x in ga]
       rb := [[(if y > 0 then y + offset else y - offset) for y in z]
                  for z in ra]
       groupPresentation(gb, rb)

   -- take quotient by removing generators specified by remgen
   quotient(a : %, remgen : List(NNI)) : % ==
       gens1 : PrimitiveArray(NNI) := a.gens
       rels1 : List(List(Integer)) := a.rels
       for toBeRemoved in remgen repeat
           gens1 := removeGen(gens1, toBeRemoved)
           rels1 := removeGen2(rels1, toBeRemoved)
       simplify(groupPresentation(entries(gens1), rels1))

   -- take quotient by adding relations specified by addrel
   quotient(a : %, addrel : List(List(Integer))) : % ==
       gens1 : List(NNI) := entries(a.gens)
       rels1 : List(List(Integer)) := a.rels
       simplify(groupPresentation(entries(gens1), concat(rels1,addrel)))

   -- directProduct of two groups
   directProduct(a : %, b : %) : % ==
       a2 : % := refactor(a)
       ga : List(NNI) := entries(a2.gens)
       ra : List(List(Integer)) := a2.rels
       sa : NNI := #ga
       b2 : % := offsetIndexes(refactor(b), sa)
       gb : List(NNI) := entries(b2.gens)
       rb : List(List(Integer)) := b2.rels
       rc : List(List(Integer)) := []
       for gax in ga repeat
           for gbx in gb repeat
               gcx : List(Integer) := [gax::Integer, gbx::Integer,
                                       gax::Integer, gbx::Integer]
               rc := concat(rc,gcx)
       groupPresentation(concat(entries(ga), gb), concat([ra, rb, rc]))

   -- cyclicGroup(n) constructs the cyclic group of order n acting
   -- on the integers 1, ..., n.
   cyclicGroup(n : PI) : % ==
       ga : List(NNI) := [1]
       ra : List(List(Integer)) := [[1 for x in 1..(n@Integer)]]
       groupPresentation(ga, ra)

   -- dihedralGroup(n) constructs the dihedral group of order 2n
   -- acting on integers 1, ..., n.
   dihedralGroup(n : PI) : % ==
       c1 : % := cyclicGroup(n)
       c2 : % := cyclicGroup(2)
       directProduct(c1, c2)

   -- symmetricGroup(n) constructs the symmetric group of order n-1.
   symmetricGroup(n : PI) : % ==
       if n<2 then return groupPresentation([])
       m : PI := (n-1) :: PI
       ga : List(NNI) := [x for x in 1..m]
       r : List(List(Integer)) := []$List(List(Integer))
       for a in 1..m repeat
           for b in 1..m repeat
               if a = b then
                   -- all generators square to 1
                   y : List(Integer) := [a,a]
                   r := cons(y, r)
               if a + 1 < b then
                   -- non squares commute
                   y : List(Integer) := [a,b,-a,-b]
                   r := cons(y, r)
               if a+1 = b then
                   -- swapping the ith and (i + 1)th position cubed
                   y : List(Integer) := [a,b,a,b,a,b]
                   r := cons(y, r)
       groupPresentation(ga, reverse!(r))

   A2D ==> TwoDimensionalArray NNI
   A1D ==> OneDimensionalArray NNI

   -- TC_state (Todd-Cox State)
   -- This holds tables while they are being built up during the
   -- Todd-Cox algorithm. Initially the elements of coset tables
   -- are set to zero then index values are incrementally
   -- derived during the algorithm. When complete the coset table
   -- will be used to determine the coset permutations.
   -- The components of this state are:
   -- coset_table          : A 2D array with
   --                      a row for each point
   --                      only 'number_of_indices' are used
   --                      other rows are empty.
   --                      a column for a generator/coset
   -- equiv_table          : Tracks coincidences
   --                      non-coincidences are self loops
   --                      coincidences are chained to non-coincidences
   -- inverse_table        : Inverse of generators (not yet used).
   -- number_of_generators : Actual number of distinct cosets.
   --                        So we can compress tables if there are
   --                        too many coincidences.
   -- number_of_indices    : allocated rows (may not be used) of
   --                        coset table.
   -- number_of_points     : size of set being permuted.
   TC_state ==> Record(coset_table : A2D, equiv_table : A1D,
                       inverse_table : A1D, closed_point : NNI,
                       number_of_generators : NNI,
                       number_of_indices : NNI, number_of_points : NNI,
                       max_number_of_indices : NNI)

   -- Local function to print CosetTable
   -- Using this avoids displaying empty rows that are not used yet.
   outCosetTable(ct : A2D,np : NNI) : OutputForm ==
       if np < 1 then return ct::OutputForm
       if nrows(ct) > np then
           part1 : NNI := np
           part2 : NNI := subtractIfCan(nrows(ct),np) :: NNI
           cts : List(A2D) := vertSplit(ct,[part1,part2])
           ct := cts.1
           --return ct::OutputForm << cts.2
       ct::OutputForm

   -- Local function to print status
   -- Using this avoids displaying empty rows that are not used yet.
   outStatus(state : TC_state) : OutputForm ==
       ct : A2D := state.coset_table
       np : NNI := state.number_of_indices
       outCosetTable(ct,np)

   -- Local function to find an index in equiv_table.
   -- Parameter information:
   -- et = equiv_table
   -- ind = index to be found
   find(et : A1D, ind : NNI) : NNI ==
       ind = 0 => ind
       qelt(et, ind) = ind => ind
       j := ind
       -- starting at ind follow pointer until j=pj
       pj := 0
       while not(j = pj) repeat
           pj := j
           j := qelt(et, j)
       j := ind
       -- track compression
       while not(j = pj) repeat
           nj := qelt(et, j)
           qsetelt!(et, j, pj)
           j := nj
       pj

   -- Local function to infer coincidencies
   -- Parameter information:
   -- ct = coset_table
   -- et = equiv_table
   -- pb = pointer begin
   -- pe = pointer end
   infer_coincidencies(ct : A2D, et : A1D, pb : NNI, pe : NNI,_
                       trace : Boolean) : Void ==
       if trace then print(message("coincidence: begin=") << pb _
                        << message(" end=") << pe)
       if pe < pb then
           (pb, pe) := (pe, pb)
       qsetelt!(et, pe, pb)
       nn := maxColIndex(ct)
       for i in 1..nn repeat
           i1 := qelt(ct, pb, i)
           i2 := qelt(ct, pe, i)
           i2 = 0 => "skip"
           if i1 = 0 then
               qsetelt!(ct, pb, i, i2)
       for i in 1..nn repeat
           i1 := qelt(ct, pb, i)
           i1 = 0 => "skip"
           i2 := qelt(ct, pe, i)
           i2 = 0 => "skip"
           i1 := find(et, i1)
           i2 := find(et, i2)
           i1 = i2 => "skip"
           infer_coincidencies(ct, et, i1, i2, trace)

   add_point ==>
        if n_ind >= o_size then
            n_ind = state.max_number_of_indices => return true
            n_size : NNI :=
                qcoerce(min(2*o_size, state.max_number_of_indices))
            n_ct_cols := maxColIndex(ct)
            nct := new(n_size, qcoerce(n_ct_cols), 0)$A2D
            net := new(n_size, 0)$A1D
            for i in 1..o_size repeat
                qsetelt!(net, i, qelt(et, i))
                for j in 1..n_ct_cols repeat
                    qsetelt!(nct, i, j, qelt(ct, i, j))
            for i in (o_size + 1)..n_size repeat
                qsetelt!(net, i, qcoerce(i))
            ct := nct
            et := net
            o_size := n_size
            state.coset_table := ct
            state.equiv_table := et
        n_ind := n_ind + 1
        if trace then
            print(message("adding action of ") << add_gen <<
                  message(" on ") << add_to << message(" to be ")
                  << n_ind)
        qsetelt!(ct, add_to, add_gen, n_ind)
        add_gen := inv_tab(add_gen)
        qsetelt!(ct, n_ind, add_gen, add_to)
        if trace then print outStatus(state)
        state.number_of_indices := n_ind


   -- Local function used by toPermutationIfCan.
   -- If we still have gaps (zeroes) in permutations then we may
   -- be able to infer a step from the relations.
   -- If there is a single step gap then we can use this to fill in
   -- information because relations are loops and therefore start and
   -- end with the same value.
   -- If permutations are changed then we return true.
   -- If so permutations is mutated.
   -- gensR is called as a Reference in case we need to
   -- eliminate duplicate points.
   inferFromRelations(state : TC_state, rels : List(List(NNI)),
                      rrels : List(List(NNI)), rel_lens : List(NNI),
                      trace : Boolean) : Boolean ==
       ct := state.coset_table
       et := state.equiv_table
       n_ind : NNI := state.number_of_indices
       inv_tab := state.inverse_table
       if trace then print(message "inferFromRelations rels=" << rels)
       add_to : NNI := 0
       add_gen : Integer := 0
       pn := state.closed_point + 1
       pn > n_ind => false
       state.closed_point := pn
       not(qelt(et, pn) = pn) => true
       o_size := maxRowIndex(ct)
       closed := false
       while not(closed) repeat
           closed := true
           not(qelt(et, pn) = pn) => break
           for rel in rels for rrel in rrels for r_len in rel_lens repeat
               -- in rel look for gap of length 1
               gb : NNI := 0 -- generator index in forward sequence
               pb : NNI := pn -- point index in forward sequence
               -- how far can we get from beginning
               i : NNI := pn
               genIn : Integer := 0
               for genIndex in rel repeat
                   genIn := genIndex
                   i := qelt(ct, i, genIndex)
                   i = 0 => break
                   i := find(et, i)
                   gb := gb + 1
                   pb := i
               gap := r_len - gb
               ge : NNI := 0 -- generator index in backward sequence
               pe : NNI := pn -- point index in backward sequence
               -- how far can we get from end
               i : NNI := pn
               genInv : Integer := 0
               for genIndex in rrel while ge < gap repeat
                   genInv := genIndex
                   i := qelt(ct, i, genIndex)
                   i = 0 => break
                   i := find(et, i)
                   ge := ge + 1
                   pe := i
               gap := gap - ge
               if gap > 1 then
                   add_to := pb
                   add_gen := genIn
                   add_point
                   gap := gap - 1
                   prel := rel
                   for i in 0..gb repeat
                       prel := rest(prel)
                   gb := gb + 1
                   while gap > 1 repeat
                       add_to := n_ind
                       add_gen := first(prel)
                       prel := rest(prel)
                       add_point
                       gap := gap - 1
                       gb := gb + 1
                   pb := n_ind
                   genIn := first(prel)
               gap = 1 =>
                   not((npe := qelt(ct, pb, genIn)) = 0) =>
                       npe = pe => "skip"
                       if trace then
                           print(message("coincidence: ") << pe <<
                                message(" ") << npe <<
                                message(" rel = ") << rel <<
                                message(" pn = ") << pn)
                       infer_coincidencies(ct, et, pe, npe, trace)
                   if trace then
                       print(message "inferFromRelations genIn=" << genIn << _
                             message " gb=" << gb)
                   not(qelt(inv_tab, genIn) = genInv) => error "impossible 3"
                   qsetelt!(ct, pb, genIn, pe)
                   not((npb := qelt(ct, pe, genInv)) = 0) =>
                       npb = pb => "skip"
                       if trace then
                          print(message("coincidence: ") << npb <<
                                message(" ") << pb <<
                                message(" rel = ") << rel <<
                                message(" pn = ") << pn)
                       infer_coincidencies(ct, et, pb, npb, trace)
                   qsetelt!(ct, pe, genInv, pb)
               gap = 0 =>
                   pb = pe => "skip"
                   -- coincidence
                   if trace then
                      print(message("coincidence: ") << pb << message(" ") <<
                          pe << message(" rel = ") << rel << message(" pn = ")
                          << pn)
                      -- print(ct::OutputForm)
                   infer_coincidencies(ct, et, pb, pe, trace)
               error "impossible 4"

           if not(closed) then
               error "impossible 5"
       true

   -- Local function used by relatorTables.
   -- Invert a map that is implemented as a 1D table
   invertMap(a : TwoDimensionalArray NNI) : TwoDimensionalArray NNI ==
       invm : TwoDimensionalArray NNI := new(nrows(a),ncols(a),0)
       for x in 1..nrows(a) repeat
           i : NNI := elt(a,x,1)
           if i ~= 0 then
               setelt!(invm,i,1,x)
       --print(message "invertMap a=" << a << message " invm=" << invm)
       invm

   -- Local function to construct relator tables.
   -- Although relator tables are used when manually calculating
   -- Todd-Cod for computations it seems easier to construct
   -- dynamically as required. So this function only used for trace
   -- output.
   relatorTables(state : TC_state,_
                 rels : List(List(Integer))) : List A2D ==
       ct := state.coset_table
       np : NNI := state.number_of_indices
       --generators : A2D := deref(gensR)
       genLists : List A2D := _
                  horizSplit(ct,ncols(ct) :: PositiveInteger)
       invGenLists : List A2D :=
                            [invertMap(a) for a in genLists]
       relators : List A2D := []
       for r in rels repeat
           relator : A2D := new(0,0,0)
           fst : Boolean := true
           for g in r repeat
               gNum : NNI := abs(g) :: NNI
               relatorn : A2D := _
                   if g>0 then genLists.gNum else invGenLists.gNum
               --print(message "relatorTables genLists=" << genLists << _
               --      message " gNum=" << gNum)
               if fst
                   then
                       relator := relatorn
                       fst := false
                   else
                       relator := horizConcat(relator,relatorn)
           relatorTrim := relator
           if nrows(relatorTrim) > np then
               part1 : NNI := np
               part2 : NNI := subtractIfCan(nrows(relatorTrim),np) :: NNI
               cts : List(A2D) := vertSplit(relatorTrim,[part1,part2])
               relatorTrim := cts.1
           relators := concat(relators,relatorTrim)
       relators

   -- convert list of generators to PermutationGroup
   generators2Permutation(state : TC_state,
                          trace : Boolean) : PermutationGroup Integer ==
       ct := state.coset_table
       if trace then print(message "generators2Permutation generators=")
       n_gens := state.number_of_generators
       n_inds := state.number_of_indices
       et := state.equiv_table
       net := new(n_inds, 0)$A1D
       j : SingleInteger := 0
       for i in 1..n_inds repeat
           not(qelt(et, i) = i) => "skip"
           j := j + 1
           qsetelt!(net, i, qcoerce(convert(j)))
       perm_lists := new(n_gens, empty())$Vector(List(Integer))
       for i in 1..n_inds repeat
           not(qelt(et, i) = i) => "skip"
           for j in 1..n_gens repeat
               kk := qelt(ct, i, j)
               kk = 0 =>
                   print(message("i = ") << i << message(" j = ") <<
                         j << message(" ct(i) = ") << row(ct, i))
                   error "incomplete coset table"
               kk := qelt(net, find(et, kk))
               qsetelt!(perm_lists, j, cons(kk, qelt(perm_lists, j)))
       pl : List(Permutation(Integer)) := []
       if trace then print(perm_lists::OutputForm)
       for j in 1..n_gens repeat
           gl := reverse!(perm_lists(j))
           p := coerceImages(gl)$Permutation(Integer)
           pl := cons(p, pl)
       pl := reverse!(pl)
       permutationGroup(pl)

   -- Convert to permutation group. Return "failed" for infinite groups.
   toPermutationIfCan(a : %) : Union(PermutationGroup Integer, "failed") ==
       toPermutationIfCan(a,false)

   -- Convert to permutation group. Return "failed" for infinite groups.
   -- This function implements the Todd-Coxeter algorithm.
   -- For more information about the algorithm see:
   -- \url{http://www.euclideanspace.com/prog/scratchpad/mycode/discrete/finiteGroup/pres2perm/}

   toPermutationIfCan(a : %,trace : Boolean
                     ) : Union(PermutationGroup Integer, "failed") ==
       toPermutationIfCan(a, [], trace)

   convert_words(words : List(List(Integer)), inv_tab : A1D
                ) : List(List(List(NNI))) ==
       nwords : List(List(NNI)) := []
       nrwords : List(List(NNI)) := []
       for word in words repeat
           nword : List(NNI) := []
           for i in word repeat
               gen :=
                   i > 0 => i
                   inv_tab(-i)
               nword := cons(qcoerce(gen)@NNI, nword)
           nwords := cons(nword, nwords)
           nrwords := cons(reverse!([inv_tab(i) for i in nword]), nrwords)
       nwords := reverse!(nwords)
       nrwords := reverse!(nrwords)
       [nwords, nrwords]

   toPermutationIfCan(a : %, sg : List(List(Integer)), trace : Boolean
                     ) : Union(PermutationGroup Integer, "failed") ==
       gs : List(NNI) := entries(a.gens)
       rs : List(List(Integer)) := a.rels
       if #gs = 0 and #rs = 0 then
           -- if no generators return trivial group
           unit : Permutation(Integer) := 1
           return permutationGroup([unit])
       if #gs > #rs then
           -- if more generators than relations then must be infinite.
           return "failed"
       numGens : NNI := #gs
       ct := new(10, 2*numGens, 0)$A2D
       et := new(10, 0)$A1D
       for i in 1..10 repeat
           qsetelt!(et, i, i)
       inv_tab := new(2*numGens, 0)$A1D
       for i in 1..numGens repeat
           qsetelt!(inv_tab, i, i + numGens)
           qsetelt!(inv_tab, i + numGens, i)
       (nrels, nrrels) := convert_words(rs, inv_tab)
       rel_lens : List(NNI) := [#rel for rel in nrels]
       loopLimit : NNI := 5000000 quo numGens
       state := [ct, et, inv_tab, 0, numGens, 1, 1, loopLimit]$TC_state
       (nsgens, nrsgens) := convert_words(sg, inv_tab)
       sgens_lens := [#word for word in nsgens]
       inferFromRelations(state, concat(nsgens, nrels),
                                   concat(nrsgens, nrrels),
                                   concat(sgens_lens, rel_lens),
                                   trace)
       -- Permutations will hold a permutation for each generator. This will
       -- be built up as we get more information.
       while true repeat
           changedByDeduction : Boolean :=
               inferFromRelations(state, nrels, nrrels, rel_lens, trace)
           if not(changedByDeduction) then
               if trace then
                   print(message("finished using ") << state.number_of_indices)
               return generators2Permutation(state, trace)
           if trace then print(message "relatorTables=" <<
                               relatorTables(state,rs))
           state.number_of_indices >= loopLimit => break
       "failed"

   -- output
   coerce(s : %) : OutputForm ==
       ps : List(NNI) := parts(s.gens)
       g : OutputForm := outputGenList(ps)
       rs : List(List(Integer)) := s.rels
       r : OutputForm := outputRelList(rs)
       hconcat([message("<"), g, message(" | "), r, message(">")])


)abbrev package GROUPPF1 GroupPresentationFunctions1
++ Description: This package implements convertion between group
++  presentation given in terms of elements of \spadtype{FreeGroup}
++  and \spadtype{GroupPresentation}.
GroupPresentationFunctions1(S : SetCategory) : with
    convert : (List(S), List(FreeGroup(S))) -> GroupPresentation
      ++ convert(lg, lr) builds group presentation from list
      ++ of generators \spad{lg} and list of relations \spad{lr}.
  == add

    convert(lg, lr) ==
        n := #lg
        nlr : List(List(Integer)) := []
        for r in lr repeat
            nr : List(Integer) := []
            fr := factors(r)
            for t in fr repeat
                k := position(t.gen, lg)
                k < 1 => error "convert: relation contains generator"
                               "not in list of generators"
                m := t.exp
                k :=
                    m < 0 => -k
                    k
                m := abs(m)
                for l in 1..m repeat
                    nr := cons(k, nr)
            nlr := cons(reverse!(nr), nlr)
        groupPresentation([l for l in 1..n], reverse!(nlr))

)if false
\section{Homology}

Intended to hold homology which is calculated using integer
Smith normal form. This is an abelian group.

It would be good if this could be modified to be based on
FreeModule by Manuel Bronstein.
)endif

)abbrev domain HOMOL Homology
++ Author: Martin Baker
++ Description:
++   Intended to hold homology which is calculated using SmithNormalForm:
++   http://www.euclideanspace.com/prog/scratchpad/mycode/topology/homology/
++ Date Created: June 2016
++ Basic Operations:
++ Related packages:
++ Related categories:
++ Related Domains: FreeModule, FiniteSimplicialComplex
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:

Homology() : Exports == Impl where
  NNI ==> NonNegativeInteger
  GENI ==> Record(vec : Vector(Integer), ord : Integer)
  SMNI ==> SmithNormalForm(Integer, _
                       Vector(Integer), _
                       Vector(Integer), _
                       Matrix(Integer))
  SRESI ==> Record(Smith : Matrix(Integer), _
                       leftEqMat : Matrix(Integer), _
                       rightEqMat : Matrix(Integer))
  x<<y ==> hconcat(x::OutputForm, y::OutputForm)

  Exports ==> SetCategory() with
    homologyGroup : (AInt : Matrix(Integer), BInt : Matrix(Integer)) -> %
      ++ construct from differential over integers
      ++ uses method described by Waldek Hebisch here:
      ++ https://groups.google.com/forum/?hl=en#!topic/fricas-devel/mLOdQ-fwbO0
    homology : (torsionVec : List(List(Integer)), torsionOrd : List(Integer),
                free1 : List(List(Integer))) -> %
      ++ construct from lists
    homology0 : () -> %
      ++ construct empty homology, useful in validation code
    homologyz : () -> %
      ++ construct Z homology, useful in validation code
    homologyzz : () -> %
      ++ construct Z*Z homology, useful in validation code
    homologyc2 : () -> %
      ++ construct C2 homology, useful in validation code
    homologyzc2 : () -> %
      ++ construct Z+C2 homology, useful in validation code
    dispGenerators : (s : %) -> OutputForm
      ++ more detailed output with generators

  Impl ==> add
   -- Representation holds torsion as vector+order
   Rep := Record(torsionPart : List(GENI), freePart : List(Vector(Integer)))

   -- construct from differential over integers
   -- AInt is input delta as matrix
   -- BInt is output delta as matrix
   -- where BInt*AInt = 0
   homologyGroup(AInt : Matrix(Integer), BInt : Matrix(Integer)) : % ==
       --
       -- validate input
       --
       --if not empty?(BInt) then
       if nrows(AInt) ~= ncols(BInt) then
           print(message("homologyGroup validation error - A rows : ") <<_
                        nrows(AInt) << message("~= B cols : ") << ncols(BInt))
       else
           zero : Matrix(Integer) := zero(nrows(BInt), ncols(AInt))
           if BInt*AInt ~= zero then
               print(message("homologyGroup validation error - B*A ~= 0 : ")
                     << BInt*AInt << message("  ~= 0 : ") << zero)
       --
       -- calculate torsion part
       --
       res : List(GENI) := empty()$List(GENI)
       smit : SRESI := completeSmith(AInt)$SMNI
       left : Matrix(Integer) := smit.leftEqMat
       m : Matrix(Integer) := smit.Smith
       leftNRows : NNI := nrows(left)
       mNRows : NNI := nrows(m)
       mNCols : NNI := ncols(m)
       for nr in 1..leftNRows repeat
           r : Vector(Integer) := row(left, nr)
           order : Integer := 1::Integer
           if nr <= mNRows and nr <= mNCols then
               order := elt(m, nr, nr)
           order <= 1 => iterate
           g : GENI := [r, order]
           res := concat(res, g)
       --
       -- calculate free part
       --
       augmented : Matrix(Integer) := vertConcat(transpose(AInt), BInt)
       --print("homologyGroup free: augmented=" << augmented)
       smitFree : SRESI := completeSmith(augmented)$SMNI
       mFree : Matrix(Integer) := smitFree.Smith
       n_rows := nrows(mFree)
       n_cols := ncols(mFree)
       kernelFree := []$List(Vector(Integer))
       for i in 1..n_cols repeat
           if i > n_rows or mFree(i, i) = 0 then
               v := new(n_cols, 0)$Vector(Integer)
               v(i) := 1
               kernelFree := cons(smitFree.rightEqMat*v, kernelFree)
       [res, reverse!(kernelFree)]

   -- construct from lists
   homology(torsionVec : List(List(Integer)), torsionOrd : List(Integer),
            free1 : List(List(Integer))) : % ==
       if #torsionVec ~= #torsionOrd then
           error "attempt to construct homology with #torsionVec ~= #torsionOrd"
       res : List(GENI) := empty()$List(GENI)
       for r1 in torsionVec for r2 in torsionOrd repeat
           r3 : GENI := [vector(r1), r2]
           res := concat(res, r3)
       kernelFree : List(Vector(Integer)) := [vector(v) for v in free1]
       [res, kernelFree]

   -- construct empty homology
   homology0() : % ==
       homology([], [], [])

   -- construct Z homology
   homologyz() : % ==
       homology([], [], [[1]])

   -- construct ZZ homology
   homologyzz() : % ==
       homology([], [], [[1, 0], [0, 1]])

   -- construct C2 homology
   homologyc2() : % ==
       homology([[1]], [2], [])

   -- construct Z+C2 homology
   homologyzc2() : % ==
       homology([[1, 0]], [2], [[0, 1]])

   -- more detailed output with generators
   dispGenerators(s : %) : OutputForm ==
       res : OutputForm := empty()$OutputForm
       s1 := s::Rep
       for g in s1.torsionPart repeat
           ln := hconcat([message("gen="), (g.vec)::OutputForm, _
               message(" ord="), (g.ord)::OutputForm])$OutputForm
           res := vconcat(res, ln)
       ln2 := hconcat([message(" free part="),
                       (s1.freePart)::OutputForm])$OutputForm
       res := vconcat(res, ln2)
       res

   -- equal if same Betti numbers and torsion coefficient
   -- This form of equality is useful for validating code. We want to check
   -- that generated homology is essentially the same as we are expecting.
   _=(a : %, b : %) : Boolean ==
       --print("homologyGroup torsionPart a : " << a.torsionPart <<_
       --      " torsionPart b : " << b.torsionPart << _
       --      " freePart a : " << a.freePart << _
       --      " numfree a : " << #(a.freePart) << _
       --      " freePart b : " << b.freePart<< _
       --      " numfree b : " << #(b.freePart))
       tora : List(GENI) := a.torsionPart
       torb : List(GENI) := b.torsionPart
       noTorsionA : Boolean := true
       noTorsionB : Boolean := true
       for ta in tora repeat
           if (ta.ord ~= 0) and (ta.ord ~= 1) then noTorsionA := false
       for tb in torb repeat
           if (tb.ord ~= 0) and (tb.ord ~= 1) then noTorsionB := false
       if noTorsionA ~= noTorsionA then return false
       #(a.freePart) = #(b.freePart)

   -- output in terms of Z (free) and C (cycles)
   -- TODO perhaps this should check if vectors are independent
   coerce(s : %) : OutputForm ==
       res : OutputForm := empty()$OutputForm
       firstTermRead : Boolean := false
       s1 := s::Rep
       --print("homologyGroup torsionPart : " << s1.torsionPart <<_
       --      " freePart : " << s1.freePart)
       nFree := #(s1.freePart)
       if nFree > 0 then
           -- TODO should check for empty list here
           res := hconcat(res, message("Z"))
           if nFree > 1 then
               res := hconcat([res, message("*"),
                               nFree::OutputForm])$OutputForm
           firstTermRead := true
       for t in s1.torsionPart repeat
           if not (t.ord = 0 or t.ord = 1) then
               if firstTermRead then
                   res := hconcat(res, message("+"))
               ln2 := hconcat([message("C"),
                               (t.ord)::OutputForm])$OutputForm
               res := hconcat(res, ln2)
               firstTermRead := true
       if not firstTermRead then
           res := hconcat(res, message("0"))
       res

--Copyright (c) 2016, Martin J Baker.
--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 Martin J Baker. 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.

)if false
\eject
\begin{thebibliography}{99}
For more details see:
[1] I have put a fuller explanation of this code here:
\url{http://www.euclideanspace.com/prog/scratchpad/mycode/discrete/finiteGroup/presentation/}

[2] Wikipedia
\url{http://https://en.wikipedia.org/wiki/Simplicial_complex/}

[3] Finite simplicial complexes in Sage
\url{http://doc.sagemath.org/html/en/reference/homology/sage/homology/simplicial_complex.html}

[4] Finite simplicial complexes in NPM
\url{https://www.npmjs.com/package/simplicial-complex}

[5] Simpcomp - a GAP package for working with simplicial complexes
\url{https://code.google.com/p/simpcomp/}

[6] A Macaulay2 package for working with simplicial complexes
\url{http://www.math.uiuc.edu/Macaulay2/doc/Macaulay2-1.8.2/share/doc/Macaulay2/SimplicialComplexes/html}

[7] Homology group method described by Waldek Hebisch here:
\url{https://groups.google.com/forum/?hl=en#!topic/fricas-devel/mLOdQ-fwbO0}

\end{thebibliography}
\end{document}
)endif