with(coxeter):
read(f_maple):

# to get the fake Coxeter number by the three 
# methods, do:
# read(fake_maple):
# then
# cox_euler(R);
# cox_sym(R);
# or
# cox_recip(R);
# where R is the name of a finite Coxeter group legal
# for the coxeter package, or one of the following,
# some of which are just ad hoc names to expand the repertoire
# past the affine groups:

# tildeA6 or similar
# tildeB6 or similar
# tildeC6 or similar
# tildeD6 or similar
# tildeE6 or similar for all other exceptional affine groups.
# ring4e5 means a graph which is a cycle with 5 elements,
# one label 4, others 3.  (Max label here = 9)
# k4_334622 means a complete graph on four vertices
# with  labels a12=3, a13=3, a14=4, a23=6, a24=2, a34=2.
# Although labels of 2 are allowed, the graph must be
# connected for this to make sense.  Max label is 9
# tildeI7 is a path with 3 vertices, one label 7 and the other 3

# To get the fake exponents, do 
#  fake_exp(R,h_proc)
# where h_proc is one of cox_euler, cox_sym or cox_recip.

# Every procedure in this file is offered with the following caveat:
# They only work for irreducible 
# Coxeter groups such that all proper parabolics are 
# either finite or irreducible or both.
# To generalize, one would need to implement 
# irreducibility-testing without relying on the coxeter 
# package's "base" and "name_of" functions,
# which only work for finite groups. 

# calculates the f_poly (reversed) but leaves 
# the Coxeter number h as a variable
# for an irreducible group of finite type or
# my infinite type.
# If this is applied to a group with a non-finite parabolic, the parabolic is 
# required to be irreducible. 
# This requires the second parameter h_proc,
# a procedure for determining a coxeter number.
# h_proc could be one of the procedures below
f_poly_h_rev:=proc(R,h_proc) local partial_f_poly,f,cm,rk; option remember;
	cm:=inf_cox_matrix(R);
	rk:=inf_rank(R);
	partial_f_poly:=0;
	f:=proc(list) ;
		try
			partial_f_poly:=partial_f_poly+f_poly_rev(name_of(linalg[submatrix](cm,list,list)));
		catch "not a finite Coxeter group":
			partial_f_poly:=partial_f_poly+f_poly_rev_inf(linalg[submatrix](cm,list,list),h_proc);
		end try;
	end proc; 
	f_on_k_sublists([seq(i,i=1..rk)],rk-1,f);
	return((h*m+2)*subs(xprime=x,int(partial_f_poly,x=0..xprime,'continuous'))/2+1);
end;

# calculates the f_poly but leaves h as a variable
# for an irreducible group of finite type or
# my infinite type.
f_poly_h:=proc(R,h_proc) option remember;
	return(x^inf_rank(R)*subs(x=1/x,f_poly_h_rev(R,h_proc)));
end;

# The euler characteristic method
#  for determining h
cox_euler:=proc(R) 
#	print(expand(subs(x=-1,factor(f_poly_h(R,cox_euler)))),expand(subs(m=m-1,x=0,factor(f_poly_h(R,cox_euler)))));
	return(solve(expand(subs(x=-1,factor(f_poly_h(R,cox_euler))))=expand(subs(m=m-1,x=0,factor(f_poly_h(R,cox_euler)))),h));
end;

# works on an irreducible Coxeter group, as long as its proper parabolics are all finite
# h_proc should be a procedure which finds a fake Coxeter number.
f_poly_rev_inf:=proc(R,h_proc) local cm,rk,partial_f_poly,f; option remember;
	cm:=inf_cox_matrix(R);
	rk:=inf_rank(R);
	partial_f_poly:=0;
	f:=proc(list) ;
		partial_f_poly:=partial_f_poly+f_poly_rev(name_of(linalg[submatrix](cm,list,list)));
	end proc; 
	f_on_k_sublists([seq(i,i=1..rk)],rk-1,f);
	return((h_proc(R)*m+2)*subs(xprime=x,int(partial_f_poly,x=0..xprime,'continuous'))/2+1);
end;

# The symmetry-based method
#  for determining h
cox_sym:=proc(R) local RmG,QmG,c,f,cm,rk; #option remember;
	cm:=inf_cox_matrix(R);
	rk:=inf_rank(R);
# calculate RmG
	RmG:=0;
	f:=proc(list) ;
		try
			RmG:=RmG+cat_formula(name_of(linalg[submatrix](cm,list,list)));
		catch "not a finite Coxeter group":
			RmG:=RmG+catalan_inf(linalg[submatrix](cm,list,list),cox_sym);
		end try;
	end proc; 
	f_on_k_sublists([seq(i,i=1..rk)],rk-1,f);
	RmG:=factor(RmG/inf_rank(R));
	#print(RmG);
# now QmG
	QmG:=factor(RmG/(m+1));
	#print(roots(QmG));
# check that this is a polynomial...
	if not(type(QmG,polynom(rational,m))) then printf("%a not a polynomial.",R); end if;
# find c  
	c:=-2*coeff(QmG,m,inf_rank(R)-3)/coeff(QmG,m,inf_rank(R)-2)/(inf_rank(R)-2);
	#print(c);
	if factor(subs(m=c-m,QmG))<>(-1)^inf_rank(R)*QmG then printf("%a not symmetric.",R); end if;
	return -2/(c+1);
end;



# finds the fake coxeter number by the reciprocity method.
# should only be tried for irreducible groups.
# This procedure leaves in the m's, rather than setting them to 1.
cox_recip:=proc(R) local partial_cat,f,partial_pos_cat,f_pos,partial_cat_in_ex,f_in_ex,h,cm,rk; option remember;
	cm:=inf_cox_matrix(R);
	rk:=inf_rank(R);
	partial_cat:=0;
	f:=proc(list) ;
		try
			partial_cat:=partial_cat+cat_formula(name_of(linalg[submatrix](cm,list,list)));
		catch "not a finite Coxeter group":
			partial_cat:=partial_cat+catalan_inf(linalg[submatrix](cm,list,list),cox_recip);
		end try;
	end proc; 
	f_on_k_sublists([seq(i,i=1..rk)],rk-1,f);
	partial_pos_cat:=0;
	f_pos:=proc(list) ;
		try
			partial_pos_cat:=partial_pos_cat+pos_cat_formula(name_of(linalg[submatrix](cm,list,list)));
		catch "not a finite Coxeter group":
			partial_pos_cat:=partial_pos_cat+pos_catalan_inf(linalg[submatrix](cm,list,list),cox_recip);
		end try;
	end proc; 
	f_on_k_sublists([seq(i,i=1..rk)],rk-1,f_pos);
	partial_cat_in_ex:=1;
	f_in_ex:=proc(list);
		if (nops(list)<rk and nops(list)>0) then 
		try
			partial_cat_in_ex:=partial_cat_in_ex+pos_cat_formula(name_of(linalg[submatrix](cm,list,list)));
		catch "not a finite Coxeter group":
			partial_cat_in_ex:=partial_cat_in_ex+pos_catalan_inf(linalg[submatrix](cm,list,list),cox_recip);
		end try;
		end if;
	end proc;
	f_on_sublists([seq(i,i=1..rk)],f_in_ex);
	#print(subs(m=1,(h*(m+1)-2)/(2*rk)*partial_pos_cat+partial_cat_in_ex));
	#print(subs(m=1,(h*m+2)/(2*rk)*partial_cat));
	return solve((h*(m+1)-2)/(2*rk)*partial_pos_cat+partial_cat_in_ex=(h*m+2)/(2*rk)*partial_cat,h);
end;

# finds the fake coxeter number by the M(G) method
# should only be tried for irreducible groups.
cox_mg:=proc(R) local cm,rk,cox,em_gee,Sigma1,Sigma2,i,f; option remember;
	cm:=inf_cox_matrix(R);
	rk:=inf_rank(R);
	Sigma1:=0;
	f:= proc(list)
		Sigma1:=Sigma1+mg(linalg[submatrix](cm,list,list))
	end proc;
	f_on_k_sublists([seq(i,i=1..rk)],rk-1,f);
	Sigma2:=0;
	f:= proc(list)
		Sigma2:=Sigma2+mg(linalg[submatrix](cm,list,list))
	end proc;
	for i from 2 to rk-1 do
		f_on_k_sublists([seq(i,i=1..rk)],i,f);
	od;
	em_gee:=Sigma1*Sigma2/(rk*(rk-1)-Sigma1);
	cox:=2/rk*(em_gee+Sigma2+rk);
	return(cox);
end proc;


#  Calculates the invariant M(G), or really M(R) for a root system R
#   For R of finite type, uses the known formulas for the Catalan number,
#   divided by m then setting m=0.
#  Otherwise, uses the method given in the paper.
mg:=proc(R) local cm,rk,em_gee,Sigma1,Sigma2,i,f; option remember;
	cm:=inf_cox_matrix(R);
	rk:=inf_rank(R);
	try
		if irreduc(name_of(R)) then
			em_gee:=rk*subs(m=0,pos_cat_formula(name_of(cm))/m);
		else em_gee:=0;
		end if;
	catch "not a finite Coxeter group":
		Sigma1:=0;
		f:= proc(list)
			Sigma1:=Sigma1+mg(linalg[submatrix](cm,list,list));
			#print(mg(linalg[submatrix](cm,list,list)));
		end proc;
		f_on_k_sublists([seq(i,i=1..rk)],rk-1,f);
		Sigma2:=0;
		f:= proc(list)
			Sigma2:=Sigma2+mg(linalg[submatrix](cm,list,list));
			#print(mg(linalg[submatrix](cm,list,list)));
		end proc;
		for i from 2 to rk-1 do
			f_on_k_sublists([seq(i,i=1..rk)],i,f);
		od;
		#print(Sigma1);
		#print(Sigma2);
		em_gee:=Sigma1*Sigma2/(rk*(rk-1)-Sigma1);
	end try;
	return(em_gee);
end proc;



# computes the "fake" exponents.
# h_proc should be a procedure which finds a fake Coxeter number.
fake_exp:=proc(R,h_proc) local cat,rts,exps,r; option remember;
	cat:=catalan_inf(R,h_proc);
	rts:=[solve(cat=0,m)];
	exps:=[];
	for r in rts do
		exps:=[op(exps),-h_proc(R)*r-1];
	od;
	return(sort(exps));
end;


# works on an irreducible Coxeter group, as long as its proper parabolics are all finite
# h_proc should be a procedure which finds a fake Coxeter number.
catalan_inf:=proc(R,h_proc) local cm,rk,partial_cat,f; option remember;
	cm:=inf_cox_matrix(R);
	rk:=inf_rank(R);
	partial_cat:=0;
	f:=proc(list) ;
		try
			partial_cat:=partial_cat+cat_formula(name_of(linalg[submatrix](cm,list,list)));
		catch "not a finite Coxeter group":
			partial_cat:=partial_cat+catalan_inf(linalg[submatrix](cm,list,list),h_proc);
		end try;
	end proc; 
	f_on_k_sublists([seq(i,i=1..rk)],rk-1,f);
	return (h_proc(R)*m+2)/(2*rk)*partial_cat;
end;

# works on an irreducible Coxeter group, as long as its proper parabolics are all finite
pos_catalan_inf:=proc(R,h_proc) local cm,rk,partial_cat,f; option remember;
	cm:=inf_cox_matrix(R);
	rk:=inf_rank(R);
	partial_cat:=0;
	f:=proc(list) ;
		partial_cat:=partial_cat+pos_cat_formula(name_of(linalg[submatrix](cm,list,list)));
	end proc; 
	f_on_k_sublists([seq(i,i=1..rk)],rk-1,f);
	return (h_proc(R)*(m+1)-2)/(2*rk)*partial_cat;
end;


# f_on_sublists applies the function f to each sublist of list
f_on_sublists:=proc(list,f) local newf, i;
	if list=[] then 
		f([]); 
		return NULL;
	end if;
	newf:=proc(x)
		f([op(1,list),op(x)]);
		f(x);	
	end proc;
	f_on_sublists([seq(op(j,list),j=2..nops(list))],newf);
	return NULL;
end;

# f_on_k_sublists applies the function f to each k-sublist of list
f_on_k_sublists:=proc(list,k,f) local newf, i;
	if list=[] then 
		if k=0 then f([]); end if; 
		return NULL;
	end if;
	newf:=proc(x)
		f([op(1,list),op(x)]);
	end proc;
	f_on_k_sublists([seq(op(j,list),j=2..nops(list))],k,f);
	f_on_k_sublists([seq(op(j,list),j=2..nops(list))],k-1,newf);
	return NULL;
end;

# gives the coxeter matrix
inf_cox_matrix:=proc(R) option remember;
	if type(R, 'array') then return R; end if;
	if substring(R,1..5)= tilde then 
		return apply(cat(substring(R,6..6),tilde),convert(substring(R,7..length(R)),decimal,10));
	elif substring(R,1..4)= ring then 
		return ring(convert(substring(R,7..length(R)),decimal,10),convert(substring(R,5..5),decimal,10));
	elif substring(R,1..2)= k4 then 
		return k4(convert(substring(R,4..4),decimal,10),convert(substring(R,5..5),decimal,10),
			convert(substring(R,6..6),decimal,10),convert(substring(R,7..7),decimal,10),
			convert(substring(R,8..8),decimal,10),convert(substring(R,9..9),decimal,10));
	else return(cox_matrix(R));
	end if;		
end;

# gives the 
inf_rank:=proc(R) option remember;
	if type(R, 'array') then return linalg[trace](R); end if;
	if substring(R,1..6)=tildeI then return 3;
	elif substring(R,1..2)= k4 then return 4;
	elif substring(R,1..5)= tilde then 
		return convert(substring(R,7..length(R)),decimal,10)+1;
	elif substring(R,1..4)= ring then 
		return convert(substring(R,7..length(R)),decimal,10);
	else return(rank(R));
	end if;		
end;


#defining the affine groups as Coxeter matrices.
# recall that the rank will be one greater than the "n"
# for n>0 although n=1 won't work in the "fake" setup.
Atilde:=proc(n)
	if n<1 then error("illegal group"); end if;
	if n=1 then return array([[1,infinity],[infinity,1]]); end if;
	return array(1..n+1,1..n+1,[seq([seq(Atilde_entry(n,i,j),j=1..n+1)],i=1..n+1)]);
end;

# for n>1.
Ctilde:=proc(n) 
	if n<2 then error("illegal group"); end if;
	return array(1..n+1,1..n+1,[seq([seq(Ctilde_entry(n,i,j),j=1..n+1)],i=1..n+1)]);
end;

# for n>2.
Btilde:=proc(n) 
	if n<3 then error("illegal group"); end if;
	return array(1..n+1,1..n+1,[seq([seq(Btilde_entry(n,i,j),j=1..n+1)],i=1..n+1)]);
end;

# for n>3.
Dtilde:=proc(n) 
	if n<4 then error("illegal group"); end if;
	return array(1..n+1,1..n+1,[seq([seq(Dtilde_entry(n,i,j),j=1..n+1)],i=1..n+1)]);
end;

Etilde:=proc(n) 
	if n=6 then
		return array(1..n+1,1..n+1,
			[	
				[1,2,3,2,2,2,2],
				[2,1,2,3,2,2,3],
				[3,2,1,3,2,2,2],
				[2,3,3,1,3,2,2],
				[2,2,2,3,1,3,2],
				[2,2,2,2,3,1,2],
				[2,3,2,2,2,2,1]
			]);
	elif n=7 then
		return array(1..n+1,1..n+1,
			[	
				[1,2,3,2,2,2,2,3],
				[2,1,2,3,2,2,2,2],
				[3,2,1,3,2,2,2,2],
				[2,3,3,1,3,2,2,2],
				[2,2,2,3,1,3,2,2],
				[2,2,2,2,3,1,3,2],
				[2,2,2,2,2,3,1,2],
				[3,2,2,2,2,2,2,1]
			]);
	elif n=8 then
		return array(1..n+1,1..n+1,
			[	
				[1,2,3,2,2,2,2,2,2],
				[2,1,2,3,2,2,2,2,2],
				[3,2,1,3,2,2,2,2,2],
				[2,3,3,1,3,2,2,2,2],
				[2,2,2,3,1,3,2,2,2],
				[2,2,2,2,3,1,3,2,2],
				[2,2,2,2,2,3,1,3,2],
				[2,2,2,2,2,2,3,1,3],
				[2,2,2,2,2,2,2,3,1]
			]);
	else error("illegal group");
	end if;
end;

Ftilde:=proc(n) 
	if n=4 then
		return array(1..n+1,1..n+1,
			[	
				[1,3,2,2,2],
				[3,1,4,2,2],
				[2,4,1,3,2],
				[2,2,3,1,3],
				[2,2,2,3,1]
			]);
	else error("illegal group");
	end if;
end;

Gtilde:=proc(n) 
	if n=2 then
		return array(1..n+1,1..n+1,
			[	
				[1,6,2],
				[6,1,3],
				[2,3,1]
			]);
	else error("illegal group");
	end if;
end;


# this is NOT an affine example,
# and the n means a label, not a rank.
# This is a path with 3 vertices, one label n, the other 3.
Itilde:=proc(n) 
	return array(1..3,1..3,
			[	
				[1,n,2],
				[n,1,3],
				[2,3,1]
			]);
end;

# a non-affine example.
# Gives a cycle of size s with one label l, others 3.
ring:=proc(s,l)
	if s<3 then error("illegal group"); end if;
	return array(1..s,1..s,[seq([seq(ring_entry(s,l,i,j),j=1..s)],i=1..s)]);
end;


# the complete graph on [4] with labels 
# a12, a13, a14, a23, a24, a34
# 
k4:=proc(a12,a13,a14,a23,a24,a34)
	return array(1..4,1..4,[[1,a12,a13,a14],[a12,1,a23,a24],[a13,a23,1,a34],[a14,a24,a34,1]]);
end;


# helpers for determining entries of the Coxeter matrices
Atilde_entry:=proc(n,i,j)
	if i=j then return 1;
	elif abs(i-j)=1 then return 3;
	elif {i,j}={1,n+1} then return 3;
	else return 2;
	end if;
end;

Ctilde_entry:=proc(n,i,j)
	if i=j then return 1;
	elif {i,j}={1,2} then return 4;
	elif {i,j}={n,n+1} then return 4;
	elif abs(i-j)=1 then return 3;
	else return 2;
	end if;
end;

Btilde_entry:=proc(n,i,j)
	if i=j then return 1;
	elif {i,j}={1,2} then return 4;
	elif {i,j}={n,n+1} then return 2;
	elif {i,j}={n,n-1} then return 3;
	elif {i,j}={n-1,n+1} then return 3;
	elif abs(i-j)=1 then return 3;
	else return 2;
	end if;
end;

Dtilde_entry:=proc(n,i,j)
	if i=j then return 1;
	elif {i,j}={1,2} then return 2;
	elif {i,j}={1,3} then return 3;
	elif {i,j}={2,3} then return 3;
	elif {i,j}={n,n+1} then return 2;
	elif {i,j}={n,n-1} then return 3;
	elif {i,j}={n-1,n+1} then return 3;
	elif abs(i-j)=1 then return 3;
	else return 2;
	end if;
end;

# helpers for determining entries of the Coxeter matrices
ring_entry:=proc(s,l,i,j)
	if i=j then return 1;
	elif abs(i-j)=1 then return 3;
	elif {i,j}={1,s} then return l;
	else return 2;
	end if;
end;

# returns the catalan number from the product formula
cat_formula:=proc(R) local cat,expo,name,first; option remember;
	if rank(R)=0 then return 1; end if;
	if irreduc(name_of(R)) then
		cat:=1;
		for expo in exponents(R) do
			cat:=cat*(m*cox_number(R)+expo+1)/(expo+1);
		od;
		return cat;
	else
		name:=name_of(R);
		first:=op(1,op(2,factors(name)));
		return (cat_formula(op(1,first))^op(2,first))*cat_formula(name/op(1,first)^op(2,first));
	end if;
end;

# returns the catalan number from the product formula
pos_cat_formula:=proc(R) local cat,expo,name,first; option remember;
	if rank(R)=0 then return 1; end if;
	if irreduc(name_of(R)) then
		cat:=1;
		for expo in exponents(R) do
			cat:=cat*(m*cox_number(R)+expo-1)/(expo+1);
		od;
		return cat;
	else
		name:=name_of(R);
		first:=op(1,op(2,factors(name)));
		return (pos_cat_formula(op(1,first))^op(2,first))*pos_cat_formula(name/op(1,first)^op(2,first));
	end if;
end;