unprotect(L,M,P,`&*`, MATRIX);

`type/M` := proc( a )
   nops(a) > 0  and { a[] } = { $1..nops(a) } 
      and type( a, list(posint) );
end:

`type/P` := proc()
   type( args, list(list(posint)) )
end:

`type/L` := proc( a )
   type( a[1], [[integer,integer],[integer,integer]] ) and
      type( a[2], prime ) and
      a[1][1][1]*a[1][2][2] mod a[2] <> a[1][1][2]*a[1][2][1] mod a[2];
end:

`type/Group` := proc()
   type( args, [list,set] )
end:

__p := 3:

ChoosePrime := proc( i::prime )
   global __p;
   __p := i;
end:

L := proc( a::list(list(integer)) )
   local matrix ;
   matrix := [ [ [a[1][1] mod __p, a[1][2] mod __p], 
                 [a[2][1] mod __p, a[2][2] mod __p] ], __p ];
   if not type( matrix, L ) then
      ERROR(`not a valid, invertible 2x2 matrix`);
   else
      RETURN(matrix);
   fi
end:

TO_P := proc( a::M )
   local unused, cycles, c;
   unused := { $1..DEGREE(a) };
   cycles := [];
   while unused <> {} do
      c := ORBT( a, unused[1] );
      unused := unused minus { op(c) };
      if nops(c) > 1 then
         cycles := [ op(cycles), c ];
      fi;
   od;
   RETURN( cycles );
end:

TO_M := proc( a::P, n::posint )
   if n < max( MOVED_POINTS(a)[] ) then
      ERROR( `degree too small` );
   else
      map( i -> a &@ i, [$1..n] );
   fi;
end:

IDENTITY_Q := proc(a)
   if type(a,M) then
      evalb( a = [ $1..DEGREE(a) ]);
   elif type(a,P) then
      evalb( a = [] );
   elif type(a,L) then
      evalb( MATRIX(a) = [[1,0],[0,1]] );
   fi;
end:   

SAME_TYPE_Q := proc(a,b)
   evalb( IDENTITY(a) = IDENTITY(b) );
end:

DEGREE := proc(a)
   if type(a,M) then 
      nops(a) ;
   elif type(a,L) then
      nops( MATRIX(a) );
   elif type(a,P) then
      ERROR(`degree is undefined`);
   fi;
end:

MATRIX := proc( a::L )
   local p;
   p := FIELD_CHAR(a);
   [ [ a[1][1][1] mod p, a[1][1][2] mod p ], 
     [ a[1][2][1] mod p, a[1][2][2] mod p ] ];
end:

FIELD_CHAR := proc( a::L )
   a[2];
end:

IDENTITY := proc(a)
   if type(a,Group) then
      IDENTITY( Generators(a)[1] );
   elif type(a,M) then    
      [ $1..DEGREE(a) ];  
   elif type(a,P) then  
      [];
   elif type(a,L) then
      [ [[1,0],[0,1]], FIELD_CHAR(a) ];
   fi;
end:

Ord := proc(a)
   local i, b;
   if type( a, Group ) then 
      RETURN( nops(Elements(a)) );
   else
      i := 1;  b:= a;
      while not IDENTITY_Q(b) do
         b := a&*b;
         i := i+1;
      od;
      RETURN(i);
   fi;
end:

CYCLE_TYPE := proc(a)
   if type(a,M) then 
      CYCLE_TYPE( TO_P(a) );
   elif type(a,P) then
      sort( map( nops, a ) );
   fi;
end:

MOVED_POINTS := proc( a::P )
   { map(op,a)[] };
end:

Inverse := proc(a)
   local i,n,v, m,p,c;
   if type(a,M) then
      n := DEGREE(a);
      v := array(1..n);
      for i from 1 to n do
         v[ a &@ i ] := i;
      od;
      RETURN( [ seq( v[i], i=1..n ) ] );
   elif type(a,P) then
      RETURN( [ seq( CYCLE_FORWARD( REVERSE(i) ), i=a ) ] );
   elif type(a,L) then
      m := MATRIX(a);
      p := FIELD_CHAR(a);
      c := 1/( m[1][1]*m[2][2] - m[1][2]*m[2][1] ) mod p;
      RETURN( [ [ [ c*m[2][2]  mod p, (-1)*c*m[1][2] mod p], 
                  [ (-1)*c*m[2][1] mod p, c*m[1][1]  mod p] ], p ] );
   else
      map( Inverse, a );
   fi;
end:

`&*` := proc(a,b)
   local n, A, B, p, c, d;
   if not SAME_TYPE_Q(a,b) then
      ERROR(`type mismatch`);
   else
      if type(a,M) then
         map( x-> a[x], b );
      elif type(a,P) then
         n := max( MOVED_POINTS(a)[], MOVED_POINTS(b)[] );
         TO_P( TO_M(a,n) &* TO_M(b,n) );
      elif type(a,L) then
         A := MATRIX(a);
         B := MATRIX(b);
         p := FIELD_CHAR(a);
         [ [ [ ( A[1][1]*B[1][1] + A[1][2]*B[2][1] ) mod p,
               ( A[1][1]*B[1][2] + A[1][2]*B[2][2] ) mod p ],
             [ ( A[2][1]*B[1][1] + A[2][2]*B[2][1] ) mod p,
               ( A[2][1]*B[1][2] + A[2][2]*B[2][2] ) mod p ] ], p ] ;
       elif type(a, set) then
          {seq( seq( c&*d, c=a), d=b)};
      fi;
   fi;
end:

Conjugate := proc(a,b)
   local c;
   if SAME_TYPE_Q(a,b) then
      a &* b &* Inverse(a);
   elif type(a,P) and ( type(b,set) or type(b,list) ) then
   { seq( Conjugate(a,c), c=b ) };
   else
      ERROR(`type mismatch`);
   fi;
end:

Commutator := proc(a,b)
   if SAME_TYPE_Q(a,b) then
      a &* b &* Inverse(a) &* Inverse(b);
   else
      ERROR(`type mismatch`);
   fi;
end:

`&@` := proc(a,i)
   local cycle, pos, j;
   if type(a,M) then
      if type(i,posint) and i <= DEGREE(a) then
         RETURN( a[i] )
      else
         ERROR(`invalid action`);
      fi;
   elif type(a,P) and type(i,posint) then
      for cycle in a do
         pos := POSITION(cycle,i);
         if pos = nops(cycle) then
            RETURN( cycle[1] );
         elif pos <> 0 then
            RETURN( cycle[ pos+1 ] );
         fi;
      od;
      RETURN(i);
   elif type(a,L) and type(i,list(integer)) and nops(i)=DEGREE(a) then
      map( x -> sum( x[j]*i[j], j=1..DEGREE(a) ) mod FIELD_CHAR(a), MATRIX(a) );
   fi;
end:

ORBT := proc(a,i)
   local o, x, v;
   if type(a,L) then 
      v := map( e -> e mod FIELD_CHAR(a), i );
   else 
      v := i;
   fi;
   o := [v];
   x := a &@ v; 
   while x <> v do
      o := [ op(o), x ];
      x := a &@ o[-1];
   od;
   o;
end:

Orbit := proc(G::Group,i)
   map( `&@`, Elements(G), i );
end:

# AUXILIARY FUNCTIONS

POSITION := proc( a::list, i )
   local j;
   for j from 1 to nops(a) do
      if i=a[j] then RETURN(j) fi;
   od;
   RETURN(0);
end:
       
REVERSE := proc( a::list )
   local i;
   [ seq( a[-i], i=1..nops(a) ) ];
end:


CYCLE_FORWARD := proc( a::list )
   [ a[-1], op( a[1..-2] ) ];
end:



GENERATE_GROUP := proc( X::set )
   local E, Y, n;
   n := nops( X[1] );
   E := {};  
   Y := X;  
   while Y <> {} do
      E := E union Y;
      Y := { seq( seq( map(i->x[y[i]],[$1..n]), y=Y ), x=X ) } minus E;
   od;
   E;
end:

# This function has multiplication built in to avoid the overhead of
# checking the type safety.

EXTEND_GROUP_M := proc( G, x )
   local E, X, Y,  n;
   if member( x, Elements(G) ) then RETURN(G) fi;
   n := nops( x );
   E := Elements( G );  
   X := [ op( Generators(G) ), x ];  
   Y := {x};
   while Y <> {} do
      E := E union 
           { seq( seq( map(i->y[g[i]],[$1..n]), y=Y ), g=Elements(G) ) };
      Y := { seq( seq( map(i->g[y[i]],[$1..n]), y=Y ), g=X) } minus E;
   od;
   [ X, E ];
end:

EXTEND_GROUP := proc( G, x )
   local E, X, Y;
   if member( x, Elements(G) ) then RETURN(G) fi;
   E := Elements( G );  
   X := [ op( Generators(G) ), x ];  
   Y := {x};
   while Y <> {} do
      E := E union { seq( seq( g&*y, y=Y ), g=Elements(G) ) }
             union { seq( seq( y&*g, y=Y ), g=Elements(G) ) };
      Y := { seq( seq( g&*y, y=Y ), g=X) } union { seq( seq( y&*g, y=Y ), g=X) } minus E;
   od;
   [ X, E ];
end:

Group := proc()
   local x, n, G, H, Y;
   for x in args do
      if not SAME_TYPE_Q(x,args[1]) then
         ERROR(`mismatched generators`);
      fi;
   od;
   if type(args[1],P) then
      n := 1;
      for x in args do
         n := max( MOVED_POINTS(x)[], n );
      od;
      G := [ [TO_M(IDENTITY(args[1]),n) ], { TO_M(IDENTITY(args[1]),n) }  ];
      for x in args do
         G := EXTEND_GROUP_M( G, TO_M(x,n) );
      od;
      G := [ map( TO_P, Generators(G) ), map( TO_P, Elements(G) ) ];
   elif type(args[1],M) then
      G := [ [IDENTITY(args[1])], {IDENTITY(args[1])} ];
      for x in args do
         G := EXTEND_GROUP_M( G, x );
      od;
   elif type(args[1],L) then
      H := {IDENTITY(args[1])};
      Y := {IDENTITY(args[1])};
      while Y <> {} do
         H := H union Y;
         Y := { seq( seq( y&*g, y=Y ), g=args ) } minus H;
      od;
      G := [ [args], H];
   fi;
   if nops( Elements(G) ) = 1 then 
      RETURN( G )
   else
      RETURN( [ G[1][2..-1], G[2] ] )
   fi;
end:

Elements   := proc( G::Group ) G[2]; end:

Generators := proc( G::Group ) G[1]; end:

ELEMENT_Q := proc( G::Group, a )
   if type(a, Group) then
      ERROR(`2nd argument cannot be a group`);
   else
      SAME_TYPE_Q(G,a) and member(a,Elements(G))
   fi;
end:

SUBGROUP_Q := proc( G::Group, H::Group )
   SAME_TYPE_Q(G,H) and nops( {Generators(H)[]} minus Elements(G) ) = 0 
end:

NORMALIZES_Q := proc( G::Group, x )
   local g ;
   if SAME_TYPE_Q(G,x) and not type(x,Group) then
      for g in Generators(G) do
         if not member( Conjugate(x,g), Elements(G) ) then
            RETURN(false);
         fi;
      od;
      RETURN(true);
   else
      ERROR(`invalid arguments`);
   fi;
end:

NORMAL_Q := proc( G::Group, H::Group )
   local i;
   if not SUBGROUP_Q(G,H) then 
      RETURN(false);
   else
      for i in Generators(G) do
         if not NORMALIZES_Q( H, x ) then RETURN(false) fi;
      od;
   fi;
   RETURN(true);
end:

CENTRALIZES_Q := proc( G::Group, x )
   local g;
   for g in Generators(G) do 
      if x &* g <> g &* x then RETURN(false) fi;
   od;
   RETURN(true);
end: 

CENTRAL_Q := proc( G::Group, H::Group )
   local x;
   if not SUBGROUP_Q(G,H) then
      RETURN(false)
   else
      for x in Generators(G) do
         if not CENTRALIZES_Q(H,x) then RETURN(false) fi;
      od;
   fi;
end:

NORMALIZER := proc( G::Group, H::Group )
   if SAME_TYPE_Q(G,H) then
      PREDICATE_SUBGROUP( G, x -> NORMALIZES_Q(H,x) );
   else
      ERROR(`groups aren't same type`)
   fi;
end:
   
CENTRALIZER := proc( G::Group, H::Group )
   if SAME_TYPE_Q(G,H) then
      PREDICATE_SUBGROUP( G, x -> CENTRALIZES_Q(H,x) );
   else
      ERROR(`groups aren't same type`)
   fi;
end:

Centre := proc( G::Group ) CENTRALIZER(G,G) end:

Stabilizer := proc( G::Group, i )
   local v;
   if ( type(IDENTITY(G),M) or type(IDENTITY(G),P) ) and type(i,posint) then
      PREDICATE_SUBGROUP( G, x ->( x&@i=i) );
   elif type(IDENTITY(G),L) and type(i,list(integer)) and 
         nops(i)=DEGREE(IDENTITY(G)) then
      v := map( e -> e mod FIELD_CHAR(IDENTITY(G)), i );
      PREDICATE_SUBGROUP( G, x -> x&@v=v );
   else
      ERROR(`action is undefined`);
   fi;
end:

PREDICATE_SUBGROUP := proc( G::Group, p )
   local L, N, a;
   N := Group( IDENTITY(G) );
   L := Elements( G );
   while L <> {} do
      a := L[1];
      if p(a) then
         N := EXTEND_GROUP(N,a);
         L := L minus Elements(N);
      else
         L := L minus map( x -> a&*x, Elements(N) );
      fi;
   od;
   N;
end:

NormalClosure := proc( G::Group, H::Group )
   local new_elts, N;
   if SAME_TYPE_Q(G,H) then
      new_elts := { seq( seq(Conjugate(g,h),h=Generators(H)), g=Generators(G) ) }
                  minus Elements(H);
      N := H;
      while new_elts <> {} do
         N := EXTEND_GROUP( H, new_elts[1] );
         new_elts := new_elts minus Elements(N);
      od;
      RETURN(N)
   else
      ERROR(`type mismatch`);
   fi;   
end:

CommutatorSubgroup := proc( G::Group )
   Group( seq( seq( Commutator(x,y), x=Generators(G) ), y=Generators(G) ) );
end:

LeftCosets := proc( G::Group, H::Group )
   local L, E;
   if SUBGROUP_Q(G,H) then
      L := [ Elements(H) ]; 
      E := Elements(G) minus Elements(H);
      while E <> {} do
         L := [ op(L), map( h-> E[1] &* h, Elements(H) ) ];
         E := E minus L[-1];
      od;
      RETURN( L );
   else
      ERROR(`the 2nd argument must be a subgroup of the 1st` );
   fi;
end:

LeftCosetReps := proc( G::Group, H::Group )
   map( X-> X[1], LeftCosets(G,H) );
end:

RightCosets := proc( G::Group, H::Group )
   map( X -> map(Inverse,X), LeftCosets(G,H) );
end:

RightCosetReps := proc( G::Group, H::Group )
   map( X -> Inverse(X[1]), LeftCosets(G,H) );
end:

ConjugacyClass := proc( G::Group, a )
   local C, new_elts;
   if ELEMENT_Q(G,a) then
      C        := {};
      new_elts := {a};
      while new_elts <> {} do
         C        := C union new_elts;
         new_elts := { seq( seq( Conjugate(x,y), x=Generators(G) ), y=new_elts) }
                     minus C;
      od;
      RETURN( C );
   else
      ERROR(`2nd argument must be an element of the 1st group`);
   fi;
end:

NumberOfConjugates := proc( G::Group, H::Group )
   if SUBGROUP_Q(G,H) then
      RETURN( Ord(G) / Ord( NORMALIZER(G,H) ) );
   else
      ERROR(`2nd argument must be an element of the 1st group`);
   fi;
end:

CycleTypes := proc( s::set)
	local typs,ts,i,p,a,b,nts;
	typs := sort( map( CYCLE_TYPE, convert(s, list) ) );
	ts := convert( typs, set);
	nts := nops(ts);
	a := [];
	for i from 1 to nts do
		member( ts[i], typs, 'p'):
		a := [op(a),p-1];
	od;
	b := [seq( a[i]-a[i-1], i=2..nts )];
	matrix( 2, nts, [ [op(ts)], [ op(b),nops(typs)-a[nts] ] ] );
end:

FLTPermutation := proc( a::L )
	local x,b,y,q,i;
	b := [];
	for x from 0 to __p-1 do
		if 	evalb( a[1][2][1]*x + a[1][2][2] mod __p <> 0 ) then
			y := (a[1][1][1]*x + a[1][1][2])/(a[1][2][1]*x + a[1][2][2]) mod __p;
		else
			y := infinity;
		fi;
		b := [op(b),y];
	od;
	if a[1][2][1] <> 0 mod __p then
		b := [op(b), a[1][1][1]/a[1][2][1] mod __p];
	else
		b := [op(b), infinity];
	fi;
	b := map( x -> x+1,b);
	member(infinity,b,'q');
	b := subsop(q=__p+1,b);
	b := TO_P( b );
	b := TO_INF( b );
	[seq( map( x->x-1,op(i,b) ),i=1..nops(b))];
end:

TO_INF := proc( b::P )
   local i,q;
   for i from 1 to nops(b) do
		if member(__p+1,op(i,b),'q') then
			RETURN(subsop( i=subsop( q=infinity,op(i,b) ),b ) );
		fi;
	od;
	RETURN(b);
end: