matskalmult := proc ( l, A )

local C;

C := map ((x,y)->x*y, Matrix(A), l);

convert ( C, listlist );

end;

#
#

matadd := proc ( A, B )

local C;

C := LinearAlgebra:-Add ( Matrix(A), Matrix(B));

convert ( C, listlist );

end;

#
#


matmult := proc ( A, B)

local C;

C := LinearAlgebra:-Multiply(Matrix(A), Matrix(B));

convert ( C, listlist );

end;

#
#

matinv := proc ( A )

local C;

C := LinearAlgebra:-MatrixInverse(Matrix(A));

convert ( C, listlist );

end;

#
#

matdet := proc ( A )

local C;

C := eval(LinearAlgebra:-Determinant (Matrix(A)));

end;

#
#

matinnerAut := proc ( A, B)

matmult ( A, matmult ( B, matinv ( A )));

end;

#######################################################################

setNN := proc ( nn )

global NN, remsN;

NN := nn;

remsN := [ seq ( i, i=0..(NN-1) ) ];

end;

#
#

setNN ( 7 );

#
#


mskalmult := proc ( l, A )

local C;

C := map ((x,y)->x*y, Matrix(A), l);

C := convert ( C, listlist );

mmod ( C );

end;

#
#

madd := proc ( A, B )

local C;

C := LinearAlgebra:-Add ( Matrix(A), Matrix(B));

C := convert ( C, listlist );

mmod ( C );

end;

#
#


mmult := proc ( A, B)

local C;

C := LinearAlgebra:-Multiply(Matrix(A), Matrix(B));

C := convert ( C, listlist );

mmod (  C  );

end;

#
#

minv := proc ( A )

local C;

C := (Inverse ( Matrix(A) ) mod NN);

C := convert ( C, listlist );

mmod ( C );

end;

#
#

mdet := proc ( A )

local C;

C := (Det (Matrix(A)) mod NN);

RETURN ( C );

mmod ( C );

end;

#
#

minnerAut := proc ( A, B)

mmult ( A, mmult ( B, minv ( A )));

end;



#
#

mmod := proc ( A )

local C;

C := Matrix(A);

convert( LinearAlgebra:-Map ((x,y)-> x mod y, C, NN), listlist );

end;

#
#

matiszero := proc ( A )

local q, A1;

q := false;

A1 := Matrix ( A );

map ( proc (x) if x <> 0 then q := true; fi; end, A1 );

not q;

end;

#
#

#Idmat := Matrix ( [[1,0],[0,1]] );

Idmat := [[1,0],[0,1]];


#
#

matsubset := proc ( lis1, lis2 )

local lis;

RETURN ( {op(lis1)} subset {op(lis2)} );

lis := [op(lis1)];

lis := map ( (xx,yy)->matmember ( xx, yy), lis, lis2 );

convert ( lis, `and` );

end;

#
#

matlisequal := proc ( lis1, lis2 )

matsubset ( lis1, lis2 ) and matsubset ( lis2, lis1 );

end;

#
#
#

matmember := proc ( A, lis )

RETURN ( matsubset ( {A}, lis ) );

lis1 := [op(lis)];

lis1 := map((xx,yy)->matiszero(matadd(matskalmult(-1,xx),yy)), lis1, A);

convert(lis1, `or` ); 

end;

#
#

matunion := proc ( lis1, lis2 )

RETURN ( {op(lis1)} union {op(lis2)} );


if nops ( lis1 ) = 0 then RETURN ( lis2 ); fi;

if nops ( lis1 ) = 1 then

	if not matmember ( lis1[1], lis2 ) then RETURN ( [op(lis1), op(lis2)] ); fi;
	RETURN ( lis2 );
	
fi;

matunion ( [lis1[1]], matunion ( lis1[2..-1], lis2 ) );

end;

#
#


gensToGroup := proc ( genlis, oper )

local frontier, res, ff, gx, newel;

frontier := genlis;

res := {op(genlis)};

while nops(frontier) > 0 do
	
	ff := frontier[1];
	frontier := frontier[2..-1];
	
	print ( "nops frontier = ", nops ( frontier ) );
	
	for gx in genlis do
	
		newel := oper ( gx, ff );
		
		if matmember ( newel, res ) then next; fi;
		
		res := res union {newel};
		
		if not matmember ( newel, frontier ) then frontier := [op(frontier),newel]; fi;
		
	od;
	
od;

[op(res)];

end;


#
#


#
#
#

generateHTilde := proc ( N )

local rems, remssqrd, ls, l, lis1, lis2, a, b ,m1, m2, lis;

rems := {op(remsN)};

remssqrd := map (x->x^2 mod N, rems );

ls := rems minus remssqrd;

if nops ( ls ) = 0 then ERROR ( "generateHTilde:: Could not find non square mod N" ); fi;

l := ls[1];

lis1 := {};
lis2 := {};

for a in remsN do
	for b in remsN do
		m1 := Matrix ( [[a, b*l],[b,a]] );
		m1 := mmod ( m1 );
		if mdet ( m1 ) <> 0 then lis1 := matunion ( {m1}, lis1 ); fi;
		m2 := Matrix ( [[a, b*l],[-b,-a]] );
		m2 := mmod ( m2 );
		if mdet ( m2 ) <> 0 then lis2 := matunion ( {m2}, lis2 ); fi;
	od;
od;

lis := matunion( lis1, lis2 );

[op(lis)];

end;

#
#


generateSl2ZmodN := proc ( )

local a, b,c, d, m;

res := [];

for a in remsN do
	for b in remsN do
		for c in remsN do
			for d in remsN do

				m := [[a,b],[c,d]];
				
				if (matdet ( m ) mod NN )= 1 then res := [op(res), m]; fi;
				
			od;
		od;
	od;
od;
	
res;


end;

#
#

generateGamma0modN := proc ( )

local a, b,c, d, m;

res := [];

for a in remsN do
	for b in remsN do
			for d in remsN do

				m := [[a,b],[0,d]];
				
				if (matdet ( m ) mod NN )= 1 then res := [op(res), m]; fi;
				
			od;
	od;
od;
	
res;


end;

#
#

generateGl2ZmodN := proc ( )

local a, b,c, d, m;

res := [];

for a in remsN do
	for b in remsN do
		for c in remsN do
			for d in remsN do

				m := Matrix ( [[a,b],[c,d]] );
				
				if mdet ( m ) <> 0 then res := [op(res), m]; fi;
				
			od;
		od;
	od;
od;
	
res;
			
end;

#
#


findInnerAutom := proc ( groupels1, groupels2, candlist )

res := [];

for x in candlist do

	tstlis := map ( (yy,xx)->minnerAut ( xx, yy ), groupels1, x );
	
	if matlisequal ( tstlis, groupels2 ) then
		
		res := [op(res), x];
	
		print ( "x found = ", x ) ; break;
		
	fi;

od;

res;

end;

#
#

# S := Matrix ( [[0,-1],[1,0]] );

# T := Matrix ( [[1,1],[0,1]] );

# T1 := matinv ( T );

S := [[0,-1],[1,0]];

T := [[1,1],[0,1]];

T1 := matinv ( T );

#
#
#

generateFundDomain := proc ( quotGroup )

local patches, frontier, ff, mm, tstlis, newel;

patches := [Idmat];

frontier := patches;

while nops ( frontier ) > 0 do

	ff := frontier[1];
	frontier := frontier[2..-1];
	
	print ( " nops(frontier) = ", nops(frontier) );
	
	for mm in [T, T1, S] do
	
		newel := matmult ( ff, mm );
		
		tstlis := map ( (xx,yy)-> matmult ( xx, yy ), patches, matinv( newel ) );
		tstlis := map ( (xx,yy) -> matmember(mmod (xx), yy), tstlis, quotGroup );
		
		if has ( tstlis, true ) then
		
			next;
			
		fi;
		
		patches := [op(patches), newel];
		
		if not matmember ( newel, frontier ) then
		
			frontier := [op(frontier), newel];
			
		fi;
				
	od;

od;


res := [];

for pel in patches do

	jlis := [];

	for mm in [T, T1, S] do
	
		newel := matmult ( pel, mm );
		
		tstlis := map ( (xx,yy)-> matmult ( xx, yy), patches, matinv ( newel ) );
		
		tstlis := map ( (xx,yy)-> matmember (mmod(xx), yy), tstlis, quotGroup );
		
		tstlis := [seq ( [i, tstlis[i]], i=1..nops(tstlis) )];
		
		jpos := select ( xx->evalb(xx[2]), tstlis )[1][1];
		
		if matmember ( newel, patches ) then
			
			jpos := -jpos;
		
		fi;
		
		jlis := [op(jlis), jpos];
		
	od;
	
	res := [op(res),[pel,jlis]];
	
od;
	
res;

end;

#
#

generateFundDomain2 := proc ( funddom, quotGroup )

local newel, res, frontier;

newel := funddom[1];

res := [newel];

frontier := [newel];

while nops ( frontier ) > 0 do

	newel := frontier[1];
	frontier := frontier[2..-1];
		
	jlist := newel[2];
	
	for j in jlist do
	
		if j > 0 then next; fi;
		
		j1 := -j;
		
		newel := funddom[j1];
		
		tstlist := map ((xx,yy)-> mmult (xx[1],yy), res, minv ( newel[1] ) );
	
		tstlist := map ((xx,yy)-> matmember ( xx, yy), tstlist, quotGroup );
		
		if has ( tstlist, true ) then next; fi;
		
				
		res := [op(res), newel];
					

		if not matmember ( newel[1], map ( xx->xx[1], frontier ) ) then
		
			frontier := [op(frontier), newel];
			
		fi;
	
	od;
		
od;

res := map ( xx -> xx[1], res );

res1 := [];

projlist := [];

for fdel in funddom do

	mfdel := fdel[1];
	
	tstlist := map ((xx,yy)->mmult(xx,yy), res, minv ( mfdel ) );
	
	tstlist := map ((xx,yy)->matmember(xx,yy), tstlist, quotGroup );
	
	tstlist := [seq([tstlist[i],i],i=1..nops(tstlist))];
	
	tstlist := select ( xx->evalb(xx[1]=true), tstlist )[1];
	
	projlist := [op(projlist), tstlist[2]];
	
od;

for pel in res do

	jlis := [];

	for mm in [T, T1, S] do
	
		newel := matmult ( pel, mm );
		
		tstlist := map ((xx,yy) -> mmult (xx,yy), res, minv(newel) );
		
		tstlist := map ((xx,yy) -> matmember(xx,yy), tstlist, quotGroup );
		
		tstlist := [seq([tstlist[i],i],i=1..nops(tstlist))];
		
		tstlist := select ( xx->xx[1] = true, tstlist );
		
		j := tstlist[1][2];
		
		if matmember ( newel, res ) then j := -j; fi;
		
		jlis := [op(jlis), j];
		
	od;
	
	res1 := [op(res1), [pel, jlis]];
		
od;	
		
[res1, projlist];

end;




#
#
#


applymatmem := [];

memcount := 1;


#
#

applymat := proc ( A, z )

global memcount, applymatmem;

local A1;

A1 := Matrix ( A );
res := simplify(limit((A1[1,1]*z1+A1[1,2]) / (A1[2,1]*z1+A1[2,2]),z1=z));

memsel := select ( (xx,yy)->
	evalb(simplify(evala(Normal(xx[1])))=simplify(evala(Normal(yy)))),
			applymatmem, res );

if nops(memsel) > 0 then

	ASSERT ( nops ( memsel ) = 1 );
	
	RETURN ( memsel[1][2] );
	
fi;

applymatmem := [op(applymatmem), [res,cat ( QQQ, memcount )]];

memcount := memcount + 1;

RETURN (applymatmem[-1][2]);

end;

#
#
#


rho := exp(2*I*Pi/3);

ival := I;

infty := infinity;

#
# generateInftyRels computes the preimages of infinity
#

generateInftyRels := proc ( funddom )

rels := {};

for pp in funddom do


	rel0 := [applymat ( pp[1], infty ), applymat ( pp[1], infty )];

	rels := rels union {rel0};

	if pp[2][1] > 0 then

		rel1 := [applymat ( matmult ( pp[1], T ), infty ), 
									applymat ( funddom[pp[2][1]][1], infty )];
		
		rels := rels union {rel1};

	fi;

	if pp[2][2] > 0 then

		rel2 := [applymat ( pp[1], infty ),
									applymat ( matmult ( funddom[pp[2][2]][1], T ), infty )];

		rels := rels union {rel2};
		
	fi;

od;

print ( rels );

rels1 := {};

# compute reflexive and symmetric closure

for rr in rels do

	rels1 := rels1 union {[rr[1],rr[2]],[rr[2],rr[1]],[rr[1],rr[1]],[rr[2],rr[2]]};
	
od;

# compute transitive Closure

tclosrels1 := transClosure ( rels1 );

print ( tclosrels1 );

equivclasses := {};

for ii in [indices(tclosrels1)] do

	equivclasses := equivclasses union {tclosrels1[op(ii)]};
	
od;

tstlist := map ((xx,yy)->[applymat(xx[1],yy),xx[1]], funddom, infty );

tstlist := [seq([op(tstlist[i]),i],i=1..nops(tstlist))];

equivres := [];

for equicl in equivclasses do
	
	equirep := select ((xx,yy)->member(xx[1],yy), tstlist, equicl);
	
	# equirep := sort ( equirep, (xx,yy)->xx[3]<yy[3] );
	
	# equirep := equirep[1][3];
	
	equirep := map ( xx->xx[3], equirep );
	
	equivres := [op(equivres), equirep];
	
od;

equivres;

end;

#
# generateRhoRels computes the preimages of rho
#

generateRhoRels := proc ( funddom )

rels := {};

for pp in funddom do

	# rel0 := [applymat ( pp[1], rho ), applymat ( matmult ( pp[1], S), rho )];

	rel0 := [applymat ( pp[1], rho ), applymat ( pp[1], rho )];

	rels := rels union {rel0};

	if pp[2][1] > 0 then

		rel1 := [applymat ( matmult ( pp[1], T ), rho ), 
									applymat ( funddom[pp[2][1]][1], rho )];
		
		rels := rels union {rel1};

	fi;

	if pp[2][2] > 0 then

		rel2 := [applymat ( pp[1], rho ),
									applymat ( matmult ( funddom[pp[2][2]][1], T ), rho )];

		rels := rels union {rel2};
		
	fi;

		
	if pp[2][3] > 0 then
								
		rel3a := [applymat ( matmult ( pp[1], S ), rho ), 
									applymat ( funddom[pp[2][3]][1], rho )];
		rel3b := [applymat ( pp[1], rho ),
									applymat ( matmult ( funddom[pp[2][3]][1], S ), rho )];

		rels := rels union {rel3a} union {rel3b};

	fi;


od;

print ( rels );

rels1 := {};

# compute reflexive and symmetric closure

for rr in rels do

	rels1 := rels1 union {[rr[1],rr[2]],[rr[2],rr[1]],[rr[1],rr[1]],[rr[2],rr[2]]};
	
od;

# compute transitive Closure

tclosrels1 := transClosure ( rels1 );

print ( tclosrels1 );

equivclasses := {};

for ii in [indices(tclosrels1)] do

	equivclasses := equivclasses union {tclosrels1[op(ii)]};
	
od;

tstlist := map ((xx,yy)->[applymat(xx[1],yy),xx[1]], funddom, rho );

tstlist := [seq([op(tstlist[i]),i],i=1..nops(tstlist))];

equivres := [];

for equicl in equivclasses do
	
	equirep := select ((xx,yy)->member(xx[1],yy), tstlist, equicl);
	
	# equirep := sort ( equirep, (xx,yy)->xx[3]<yy[3] );
	
	# equirep := equirep[1][3];
	
	equirep := map ( xx->xx[3], equirep );
	
	equivres := [op(equivres), equirep];
	
od;

equivres;

end;

#
# generateIRels computes the preimages of I
#

generateIRels := proc ( funddom )

rels := {};

for pp in funddom do


	rel0 := [applymat ( pp[1], ival ), applymat ( pp[1], ival )];

	rels := rels union {rel0};

		
	if pp[2][3] > 0 then
								
		rel3a := [applymat ( matmult ( pp[1], S ), ival ), 
									applymat ( funddom[pp[2][3]][1], ival )];
		rel3b := [applymat ( pp[1], ival ),
									applymat ( matmult ( funddom[pp[2][3]][1], S ), ival )];

		rels := rels union {rel3a} union {rel3b};

	fi;


od;

print ( rels );

rels1 := {};

# compute reflexive and symmetric closure

for rr in rels do

	rels1 := rels1 union {[rr[1],rr[2]],[rr[2],rr[1]],[rr[1],rr[1]],[rr[2],rr[2]]};
	
od;

# compute transitive Closure

tclosrels1 := transClosure ( rels1 );

print ( tclosrels1 );

equivclasses := {};

for ii in [indices(tclosrels1)] do

	equivclasses := equivclasses union {tclosrels1[op(ii)]};
	
od;

tstlist := map ((xx,yy)->[applymat(xx[1],yy),xx[1]], funddom, ival );

tstlist := [seq([op(tstlist[i]),i],i=1..nops(tstlist))];

equivres := [];

for equicl in equivclasses do
	
	equirep := select ((xx,yy)->member(xx[1],yy), tstlist, equicl);
	
	# equirep := sort ( equirep, (xx,yy)->xx[3]<yy[3] );
	
	# equirep := equirep[1][3];
	
	equirep := map ( xx->xx[3], equirep );
	
	equivres := [op(equivres), equirep];
	
od;

equivres;

end;

#
#
#
#

transClosure := proc ( rels )

local elems, elemfollowers, el;

elems := {};

elems := `union` (op ( [seq ( {rels[i][1],rels[i][2]}, i=1..nops(rels) )] ));

elemfollowers := table ( [] );

for el in elems do
	
	elemfollowers[el] := {};
	
od;

for rel in rels do

	start := rel[1];
	ziel := rel[2];
	
	elemfollowers[start] := elemfollowers[start] union elemfollowers[ziel] union {ziel};
	
	for el in elems do
	
		if member ( start, elemfollowers[el] ) then
		
			elemfollowers[el] := elemfollowers[el] union elemfollowers[start];
			
		fi;
		
	od;
	
od;

eval(elemfollowers);

end;

#
#