with(coxeter):
with(weyl):

#  to check Theorem {th:restrict}  2.7, do:
#  read(compat_maple):
#  for R in [H3,H4,F4,E6,E7,E8] do check_res_all(R); od;
#  It should take about ?? hours


#gives a bipartition of the diagram
#in the form of two sets of indices
bipartition:=proc(R) local i,neighbors,j,bip;
	if rank(R)=0 then return([[],[]]); end if;
	for i from 1 to rank(R) do
		neighbors:=[];
		for j from 1 to rank(R) do
			if not(i=j) and not(fuzzy_orth(op(i,base(R)),op(j,base(R)))) then 
				neighbors:=[op(neighbors),j];
			end if; 
			if nops(neighbors)>1 then 
				break;
			end if;
		od;
		if neighbors=[] then
			bip:=bipartition(subsop(i=NULL,base(R)));
			bip:=shift_bip(bip,i);
			return([[i,op(op(1,bip))],op(2,bip)]);
		elif nops(neighbors)=1 then
			bip:=bipartition(subsop(i=NULL,base(R)));
			bip:=shift_bip(bip,i);
			if op(1,neighbors) in op(1,bip) then
				return([op(1,bip),[i,op(op(2,bip))]]);
			else
				return([[i,op(op(1,bip))],op(2,bip)]);
			end if;
		end if;
	od;
end;


#  A helper function which lets the bipartition
#  of a parabolic subgroup embed in a bipartition of 
#  the big group.
shift_bip:=proc(bip,i) local p1,p2,j;
	p1:=op(1,bip);
	p2:=op(2,bip);
	for j from 1 to nops(p1) do
		if op(j,p1)>= i then 
			p1:=subsop(j=op(j,p1)+1,p1);
		end if;
	od;
	for j from 1 to nops(p2) do
		if op(j,p2)>= i then 
			p2:=subsop(j=op(j,p2)+1,p2);
		end if;
	od;
	return([p1,p2]);
end;

# the maps tau_plus and tau_minus from Generalized Associahedra
# according to the bipartition defined in bipartition(R).
# usage tau_plus(R,alpha) or tau_plus(R,alpha,bip)
# bip is a bipartition (a list of two lists of indices of base(R)
# which bipartitions the diagram.)
tau_plus:=proc() local i,mi,R,alpha,bip;
	R:=args[1];
	alpha:=args[2];
	if nargs=3 then bip:=args[3] else bip:=bipartition(R); end if;
	mi:=op(1,bip);
	for i in mi do
		if fuzzy_eq(op(i,base(R)),(-1)*alpha) then
			return alpha;
		end if;
	od;
	return reflect(seq(op(i,base(R)),i=op(2,bip)),alpha);
end;	

tau_minus:=proc() local i,p,R,alpha,bip;
	R:=args[1];
	alpha:=args[2];
	if nargs=3 then bip:=args[3] else bip:=bipartition(R); end if;
	p:=op(2,bip);
	for i in p do
		if fuzzy_eq(op(i,base(R)),(-1)*alpha) then 
			return alpha;
		end if;
	od;
	return reflect(seq(op(i,base(R)),i=op(1,bip)),alpha);
end;	

# rotation, forward and reverse.
# same usage as tau_plus and tau_minus
rot:=proc() local R,alpha,bip;
	R:=args[1];
	alpha:=args[2];
	if nargs=3 then bip:=args[3] else bip:=bipartition(R); end if;
	return tau_plus(R,tau_minus(R,alpha,bip),bip);
end;


# col_rt is a "colored root" [alpha,i]
# takes three parameters R,m,col_rt
# with an optional 4th parameter bip, a bipartition.
m_rot:=proc() local R,m,col_rt,bip,alpha,i,b;
	R:=args[1];
	m:=args[2];
	col_rt:=args[3];
	if nargs=4 then bip:=args[4] else bip:=bipartition(R); end if;
	alpha:=op(1,col_rt);
	i:=op(2,col_rt);
	if i=m then return [rot(R,alpha,bip),1]; end if;
	for b in base(R) do
		if fuzzy_eq(b,(-1)*alpha) then 
			return [rot(R,alpha,bip),1];
		end if;
	od;
	return [alpha,i+1];
end;

# usual compatibility as defined by Fomin and Zelevinsky
compat:=proc() local i,R,alpha1,alpha2,bip; option remember;
	R:=args[1];
	alpha1:=args[2];
	alpha2:=args[3];
	if nargs=4 then bip:=args[4] else bip:=bipartition(R); end if;
	for i from 1 to rank(R) do
		if fuzzy_eq(op(i,base(R)),(-1)*alpha1) then
			return fuzzy_orth(op(i,weights(R)),alpha2);
		end if;
	od;
	return(compat(R,rot(R,alpha1,bip),rot(R,alpha2,bip),bip));
end;

# compatibility in the "m" version
m_compat:=proc() local R,m,col_rt1,col_rt2,bip,alpha1,i; option remember;
	R:=args[1];
	m:=args[2];
	col_rt1:=args[3];
	col_rt2:=args[4];
	if nargs=5 then bip:=args[5] else bip:=bipartition(R); end if;
	for i from 1 to rank(R) do
		if fuzzy_eq(op(i,base(R)),(-1)*op(1,col_rt1)) then
			return fuzzy_orth(op(i,weights(R)),op(1,col_rt2));
		end if;
	od;
	return(m_compat(R,m,m_rot(R,m,col_rt1,bip),m_rot(R,m,col_rt2,bip),bip));
end;

#  This is by the definition.
#  Previous is by the proposition.
m_compat_alt:=proc() local R,m,col_rt1,col_rt2,bip,alpha1,alpha2,i1,i2; option remember;
	R:=args[1];
	m:=args[2];
	col_rt1:=args[3];
	col_rt2:=args[4];
	if nargs=5 then bip:=args[5] else bip:=bipartition(R); end if;
	alpha1:=op(1,col_rt1);
	alpha2:=op(1,col_rt2);
	i1:=op(2,col_rt1);
	i2:=op(2,col_rt2);
	if i1>i2 and dist(R,alpha1,bip) <= dist(R,alpha2,bip) then
		return compat(R,rot(R,alpha1,bip),alpha2,bip);
	end if;
	if i1<i2 and dist(R,alpha1,bip) >= dist(R,alpha2,bip) then
		return compat(R,alpha1,rot(R,alpha2,bip),bip);
	end if;
	return compat(R,alpha1,alpha2,bip);
end;

# the distance of a root from a negative simple
# in the ordinary (non-m) rotation
# takes parameters R and alpha with an optional bip, a bipartition.
dist:=proc() local R,alpha,bip,r; option remember;
	R:=args[1];
	alpha:=args[2];
	if nargs=3 then bip:=args[3] else bip:=bipartition(R); end if;
	for r in base(R) do
		if iprod(r+alpha,r+alpha) < `coxeter/default`[epsilon] then return 0; end if;
	od;
	return dist(R,rot(R,alpha,bip))+1;
end;

# does the function f on every colored root
f_on_coloreds:=proc(R,m,f) local a,b,i;
	for a in base(R) do f([-a,1]); od;
	for b in pos_roots(R) do
		for i from 1 to m do f([b,i]); od;
	od;
	return null;
end;

list_coloreds:=proc(R,m) local l,f;
	l:=[];
	f:=proc(x) l:=[op(l),x]; end proc;
	f_on_coloreds(R,m,f);
	return(l);
end;
	
#  applies the function f to every simplex with k+1 vertices.
#  each simplex is a set of roots
f_on_k_simplices:= proc(R,m,k,f)
	f_on_k_simp_helper(R,m,{},list_coloreds(R,m),k,f);
end;

f_on_k_simp_helper:=proc(R,m,so_far,list,k,f) local i;
	if k=0 then f(so_far); return NULL; end if;
	if list=[] then return NULL; end if;
	if m_compat_list(R,m,op(1,list),so_far) then
		f_on_k_simp_helper(R,m,so_far union {op(1,list)},[seq(op(j,list),j=2..nops(list))],k-1,f);
	end if;
	f_on_k_simp_helper(R,m,so_far,[seq(op(j,list),j=2..nops(list))],k,f);
	return NULL;
end;

m_compat_list:=proc(R,m,alpha,list) local ans,beta;
	ans:=true;
	for beta in list do
		ans:=ans and m_compat(R,m,alpha,beta);
	od;
	return ans;
end;

# lists k simplices
list_simp:=proc(R,m,k) local l,f;
	l:=[];
	f:=proc(x) l:=[op(l),x]; end proc;
	f_on_k_simplices(R,m,k,f);
	return(l);
end;
	
# set of k simplices
set_simp:=proc(R,m,k) local l,f;
	l:={};
	f:=proc(x) l:=l union {x}; end proc;
	f_on_k_simplices(R,m,k,f);
	return(l);
end;

 	

# counts k-simplices
num_simp:=proc(R,m,k) local num,f;
	num:=0;
	f:=proc(x) num:=num+1; end proc;
	f_on_k_simplices(R,m,k,f);
	return num;
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;

# does this for all i and prints a report
check_res_all:=proc(R) local i;
	interface(quiet=true);
	for i from 1 to rank(R) do
		print(R,i,check_restriction(R,i));
	od;
	interface(quiet=false);
end proc;

#  Checks that two colored roots in a maximal parabolic are
#  compatible in the parabolic if and only if they are 
#  compatible in the big group.  I denotes the generator missing
#  in the parabolic.
#  Just the m=2 case, which is sufficient.
check_restriction:=proc(R,i) local RPrime,ans,bip,bipPrime,col_rts,d1,d2;
	ans:=true;
	RPrime:=subsop(i=NULL,base(R));
	bip:=bipartition(R);
	bipPrime:=collapse_bip(bip,i);
	#print(i,bip,bipPrime);
	col_rts:=list_coloreds(RPrime,2);
	for d1 in col_rts do
		for d2 in col_rts do
			ans:=ans and not(m_compat(R,2,d1,d2,bip) xor m_compat(RPrime,2,d1,d2,bipPrime));
			#print(d1,d2,m_compat(R,2,d1,d2,bip),m_compat(RPrime,2,d1,d2,bipPrime));
			if (m_compat(R,2,d1,d2,bip) xor m_compat(RPrime,2,d1,d2,bipPrime)) then
				print(d1,d2,m_compat(R,2,d1,d2,bip),m_compat(RPrime,2,d1,d2,bipPrime));
			end if;
		od;
	od;
	return(ans);
end;

# a procedure to make a bipartition for the parabolic
# obtained by deleting the ith entry of base(R).
# bip is a bipartition for R.
collapse_bip:=proc(bip,i) local j,p1,p2;
	p1:=op(1,bip);
	p2:=op(2,bip);
	for j from 1 to nops(p1) do
		if op(j,p1)=i then 
			p1:=subsop(j=NULL,p1);
			break;
		end if;
	od;
	for j from 1 to nops(p2) do
		if op(j,p2)=i then 
			p2:=subsop(j=NULL,p2);
			break;
		end if;
	od;
	for j from 1 to nops(p1) do
		if op(j,p1)>i then 
			p1:=subsop(j=op(j,p1)-1,p1);
		end if;
	od;
	for j from 1 to nops(p2) do
		if op(j,p2)>i then 
			p2:=subsop(j=op(j,p2)-1,p2);
		end if;
	od;
	return([p1,p2]);
end;



# equality of vectors up to the coxeter package's epsilon.
fuzzy_eq:=proc(a,b)
	return(evalb(iprod(a-b,a-b)<`coxeter/default`[epsilon]^2));
end;

# orthogonality of vectors up to the epsilon
fuzzy_orth:=proc(a,b)
	return(evalb(abs(iprod(a,b))<`coxeter/default`[epsilon]^2));
end;

# positivity of scalars up to epsilon
fuzzy_pos:=proc(a)
	return(evalb(a>`coxeter/default`[epsilon]));
end;

# positivity of scalars up to epsilon
fuzzy_neg:=proc(a)
	return(evalb(a<-`coxeter/default`[epsilon]));
end;

# zero-ness of scalars up to epsilon
fuzzy_zero:=proc(a)
	return(evalb(abs(a)<`coxeter/default`[epsilon]));
end;