#  Maple routines for sortable elements for finite Coxeter groups.
with(coxeter):
with(posets):
with(weyl):


#  IMPORTANT:
#  To improve speed in the crystallographic case, where 
#   everything can be done in exact arithmetic, there is a global 
#  boolean variable crys.  
#  If crys is true, then these programs run routines that are only valid
#  in the crystallographic case.
#  If it's false, these programs run routines that are valid in either case.
#  By default, crys:=true, so if you want to run one of these checks 
#  in the noncrystallographic case, you must first set crys:=false; 
crys:=true;

# checks Proposition 3.1  {family}  
# for each choice of coxeter element and returns
# the final result (using AND)
check_family_all:=proc(R) local f, num_checked,results,inttemp;
	results:=true;
	num_checked:=0;
	inttemp:=interface(quiet);
	interface(quiet=true);
	f:=proc(cox)  local restmp,t;
		t:=time();
		restmp:=check_family(R,cox);
		results:=results and restmp; 
		num_checked:=num_checked+1;
		printf("%a    %a    %a    %a       %a\n",name_of(R),num_checked,cox,restmp,time()-t);
	end proc;
	f_on_cox(R,f);
	interface(quiet=inttemp);
	return(results);
end proc;

# Checks Proposition {family} for a fixed choice of coxeter element
# as stated in the text, it's only necessary to check (ii)
# The test, as programmed is not valid in rank two, because it
# relies on the fact, (true in higher ranks) that the only labels
# are 2,3,4,5.  Thus there is a finite set of angles which characterize
# when two reflections are the canonical generators of an
#  irreducible rank-two parabolic.
#
#  once you have a pair t1,t2 of reflections canonically generating an
#  irreducible rank-two parabolic, for each initial s, the property we 
#  want to check is:
#  if neither of them equals s, then st1s and st2s should appear in the 
#  same sequence in the total order with respect to cox as in the 
#  total order with respect to s.cox.s
#  if one of them is s, then the order should be reversed. 
check_family:=proc(R,cox)  local order,len,ans,s,coxprime,orderprime,i,j,ang,refli,reflj,little_ans ;
	ans:=true;
	order:=ref_order(R,cox);
	#print(order);
	len:=nops(order);
	for s from 1 to rank(R) do
		if nops(reduce([s,op(cox)],R))<nops(cox) then
			coxprime:=reduce([s,op(cox),s],R);
			orderprime:=ref_order(R,coxprime);
			#print(s,orderprime);
			for i from 1 to len do
				for j from i+1 to len do
					ang:=iprod(op(i,order),op(j,order))
						/sqrt(iprod(op(i,order),op(i,order))*iprod(op(j,order),op(j,order)));
					# this tests if the two canonically generate an irreducible rank-two parabolic
					if crys then
						if ang=-1/2 or ang=-sqrt(2)/2 then
							refli:=reflect(op(s,base(R)),op(i,order));
							reflj:=reflect(op(s,base(R)),op(j,order));
							little_ans:=pm_in_order(refli,reflj,orderprime);
							# takes care of "reversing" the order if t1 = s
							if op(i,order)=op(s,base(R)) then little_ans:=not(little_ans); end if;
							ans:=ans and little_ans;
						end if;
					else
						if fuzzy_eq(evalf(ang),-1/2) or fuzzy_eq(evalf(ang),evalf(-sqrt(2)/2)) or
						             fuzzy_eq(evalf(ang),-0.8090169941) then
							refli:=reflect(op(s,base(R)),op(i,order));
							reflj:=reflect(op(s,base(R)),op(j,order));
							little_ans:=fuzzy_pm_in_order(refli,reflj,orderprime);
							# takes care of "reversing" the order if t1 = s
							if fuzzy_vec_eq(op(i,order),op(s,base(R))) then 
								little_ans:=not(little_ans); 
							end if;
							ans:=ans and little_ans;
						end if;
					end if;
				od;
			od;
		end if;
	od;
	return(ans);
end proc;


# checks Proposition 3.2 {restrict or} 
# for each choice of coxeter element and returns
# the final result (using AND)
check_res_all:=proc(R) local f, num_checked,results,inttemp;
	results:=true;
	num_checked:=0;
	inttemp:=interface(quiet);
	interface(quiet=true);
	f:=proc(cox)  local restmp,t;
		t:=time();
		restmp:=check_res(R,cox);
		results:=results and restmp; 
		num_checked:=num_checked+1;
		printf("%a    %a    %a    %a           %a\n",name_of(R),num_checked,cox,restmp,time()-t);
	end proc;
	f_on_cox(R,f);
	interface(quiet=inttemp);
	return(results);
end proc;

# Checks Proposition {restrict or} for a fixed choice of coxeter element
# This makes the reflection ordering on positive roots for R
#  using ref_order(R,cox), then does all possible deletions of a single 
#  generator.  For each one, it gets a reflection order and checks that 
#  no pair of noncommuting generators is ordered opposite to the 
#  ordering on the larger root system.
check_res:=proc(R,cox) local big_order,small_order,Rpara,i,j,k,ans;
	ans:=true;
	big_order:=ref_order(R,cox);
	for i from 1 to rank(R) do
		Rpara:=subsop(i=NULL,base(R));
		small_order:=ref_order(Rpara,res_cox(cox,i));
		#printf("%a\n",small_order);
		#printf("%a\n",big_order);
		if crys then
			for j from 1 to nops(small_order) do
	 			for k from j+1 to nops(small_order) do
					print(op(j,small_order),op(k,small_order));
					ans:=ans and not_misordered(op(j,small_order),op(k,small_order),big_order);
				od;
			od;
			#printf("%a   %a   %a\n",name_of(R),cox,i);
		else 
			for j from 1 to nops(small_order) do
	 			for k from j+1 to nops(small_order) do
					print(op(j,small_order),op(k,small_order));
					ans:=ans and fuzzy_not_misordered(op(j,small_order),op(k,small_order),big_order);
				od;
			od;
			#printf("%a   %a   %a\n",name_of(R),cox,i);
		end if;
	od;
	return ans;
end proc;

# checks the claim in the proof of Lemma 4.6 {align sc}
# for each choice of coxeter element and returns
# the final result
check_claim_all:=proc(R) local f, num_checked,results,inttemp;
	results:=true;
	num_checked:=0;
	inttemp:=interface(quiet);
	interface(quiet=true);
	f:=proc(cox)  local restmp,t;
		t:=time();
		restmp:=check_claim(R,cox);
		results:=results and restmp; 
		num_checked:=num_checked+1;
		printf("%a    %a    %a    %a                  %a\n",name_of(R),num_checked,cox,restmp,time()-t);
	end proc;
	f_on_cox(R,f);
	interface(quiet=inttemp);
	return(results);
end proc;

# Checks the claim in the proof of Lemma {align sc}  that says that 
#  if s is initial in c and t2 is not in W_<s> then there exists
#  t1 preceding t2 in the reflection order such that t1 is a canonical generator
#  of the rank-two parabolic they generate, while t2 is not a canonical generator,
#  and with the property that t1 is also not in W_<s>.
#  In fact, we only need to check the following:
#  there exists t1 preceding t2 with t1 not in W_<s> such that
#  reflecting the root for t1 with respect to t2 and negating
#  gives a root later in the ordering.  (If so, then t1 or some earlier root
#  is the canonical generator of the parabolic subgroup containing t1 and t2
#  and t2 is not canonical.)
check_claim:=proc(R,cox) local i,j,k,l,m,wei,ans,local_ans,rt,ord;
	wei:=weights(R);
	ans:=true;
	ord:=ref_order(R,cox);
	for j in firsts(R,cox) do
		if crys then
			# find a root (not equal to the root for s) in the ordering that is not in W_<s>
			for k from 1 to nops(ord) do
				if not(op(k,ord)=op(j,base(R))) and not(orth(op(k,ord),op(j,wei))) then
					local_ans:=false;
					#  check all roots before in the ordering
					for l from 1 to k-1 do
						#  finds one not in W_<s> which reflects t1 to the negation of something later.
						if not(orth(op(l,ord),op(j,wei))) then
							rt:=(-1)*reflect(op(k,ord),op(l,ord));
							printf("%a   %a   %a\n\n",k,l,rt);
							for m from k+1 to nops(ord) do
								if op(m,ord)=rt then 
									printf("%a   %a\n",m,rt);
									local_ans:=true; 
									break;
								end if;
							od;
						if local_ans=true then break; end if;
						end if;
					od;
					if local_ans=false then print(`lemma failed`,j,k); end if;
					ans:=ans and local_ans;
				end if;
			od;
			#printf("%a   %a   %a\n",R,cox,j);
		else
			for k from 1 to nops(ord) do
				if not(fuzzy_vec_eq(op(k,ord),op(j,base(R)))) and not(fuzzy_orth(op(k,ord),op(j,wei))) then
					local_ans:=false;
					for l from 1 to k-1 do
						if not(fuzzy_orth(op(l,ord),op(j,wei))) then
							rt:=(-1)*reflect(op(k,ord),op(l,ord));
							printf("%a   %a   %a\n\n",k,l,rt);
							for m from k+1 to nops(ord) do
								if fuzzy_vec_eq(op(m,ord),rt) then 
									printf("%a   %a\n",m,rt);
									local_ans:=true; 
									break;
								end if;
							od;
						if local_ans=true then break; end if;
						end if;
					od;
					if local_ans=false then print(`lemma failed`,j,k); end if;
					ans:=ans and local_ans;
				end if;
			od;
			#printf("%a   %a   %a\n",R,cox,j);
		end if;
	od;
	return(ans);
end proc;

# checks Lemma 6.6 {nc lemma}
# for each choice of coxeter element and returns
# the final result
check_nc_lemma_all:=proc(R) local f, results,num_checked,inttemp;
	results:={};
	num_checked:=0;
	inttemp:=interface(quiet);
	interface(quiet=true);
	f:=proc(cox)  local restmp,t;
		t:=time();
		restmp:=check_nc_lemma(R,cox);
		results:=results union {restmp}; 
		num_checked:=num_checked+1;
		printf("%a    %a    %a    %a       %a\n",name_of(R),num_checked,cox,restmp,time()-t);
	end proc;
	f_on_cox(R,f);
	interface(quiet=inttemp);
	return(results);
end proc;

#  Checks Lemma~\ref{nc lemma}
check_nc_lemma:=proc(R,cox) local firsts_lasts,ans,f,wei,B;
	B:=base(R);
	ans:=true;
	wei:=weights(R);
	firsts_lasts:=[op(firsts(R,cox)),op(lasts(R,cox))];
	# this will be a function run on all sortable elements
	# using f_on_sortas_and_base_images
	if crys then
		f:=proc(word,base_images,vec) local i,j,k,cov_known,cov;
			#print(word,base_images);
			cov_known:=false;
			for i in firsts_lasts do
				if member(-op(i,base_images),B,'k') then
					if not cov_known then
						cov:=cov_refs(R,word,vec); 
						cov_known:=true;
					end if;
					#print(i,k,cov);
					for j from 1 to nops(cov) do 
						ans:=ans and (op(i,B)=-op(j,cov) or orth(op(i,wei),op(j,cov)));
					od;
				end if;
			od;
		end proc;
	else 
		f:=proc(word,base_images,vec) local i,j,k,cov_known,cov;
			#print(word,base_images);
			cov_known:=false;
			for i in firsts_lasts do
				if fuzzy_vec_member(-op(i,base_images),B,'k') then
					if not cov_known then
						cov:=cov_refs(R,word,vec); 
						cov_known:=true;
					end if;
					#print(i,k,cov);
					for j from 1 to nops(cov) do 
						ans:=ans and (fuzzy_vec_eq(op(i,B),-op(j,cov))
						       or fuzzy_orth(op(i,wei),op(j,cov)));
					od;
				end if;
			od;
		end proc;
	end if;
	f_on_sortas_and_base_images(R,cox,f);
	return(ans);
end proc;

# checks Lemma 6.7 {nc lemma 2}
# for each choice of coxeter element and returns
# the list of all different results.
check_nc_lemma2_all:=proc(R) local f, num_checked,results,inttemp;
	results:=true;
	num_checked:=0;
	inttemp:=interface(quiet);
	interface(quiet=true);
	f:=proc(cox)  local restmp,t;
		t:=time();
		restmp:=check_nc_lemma2(R,cox);
		results:=results and restmp; 
		num_checked:=num_checked+1;
		printf("%a    %a    %a    %a         %a\n",name_of(R),num_checked,cox,restmp,time()-t);
	end proc;
	f_on_cox(R,f);
	interface(quiet=inttemp);
	return(results);
end proc;

# Checks lemma~\ref{nc lemma 2} for a fixed choice of coxeter element.
check_nc_lemma2:=proc(R,cox) local ans,B,int,la,f;
	ans:=true;
	B:=base(R);
	int:=interior_pt(R);
	la:=lasts(R,cox);
	if crys then
		f:=proc(sort,base_images) local i;
			#print(sort);
			for i in la do
				#print(op(i,base_images),iprod(int,op(i,base_images)),member(-op(i,base_images),B));
				if iprod(int,op(i,base_images))<0 then
					ans:=ans and member(-op(i,base_images),B);
				end if;
			od;
		end proc;
	else
		f:=proc(sort,base_images) local i;
			#print(sort);
			for i in la do
				#print(op(i,base_images),iprod(int,op(i,base_images)),fuzzy_vec_member(-op(i,base_images),B));
				if fuzzy_neg(iprod(int,op(i,base_images))) then
					ans:=ans and fuzzy_vec_member(-op(i,base_images),B);
				end if;
			od;
		end proc;
	end if;
	f_on_sortas_and_base_images(R,cox,f);
	return(ans);
end proc;

# checks Proposition {compat} 7.2
# for each choice of coxeter element and returns
# the list of all different results.
check_compat_all:=proc(R) local f, num_checked,results,inttemp;
	results:=true;
	num_checked:=0;
	inttemp:=interface(quiet);
	interface(quiet=true);
	f:=proc(cox)  local restmp,t;
		t:=time();
		restmp:=check_compat(R,cox);
		results:=results and restmp; 
		num_checked:=num_checked+1;
		printf("%a    %a    %a    %a         %a\n",name_of(R),num_checked,cox,restmp,time()-t);
	end proc;
	f_on_cox(R,f);
	interface(quiet=inttemp);
	return(results);
end proc;

#  This checks Proposition {compat}
#  by verifying that the c-compatibility relation
#  given by compat_rel_cl satisfies conditions (i) and (ii) of
#  Proposition {compat}
check_compat:=proc(R,cox)
	return(check_compat_i(R,cox) and check_compat_ii(R,cox));
end proc;

# checks condition (i)
check_compat_i:=proc(R,cox) local ans,i,pair,pairs,Rprime,t;
	ans:=true;
	for i from 1 to rank(R) do
		print(`compat (i)`,i);
		pairs:={};
		for pair in compat_rel_cl(R,cox) do
			if [-i] in pair then
				pairs:=pairs union {pair};
			end if;
		od;
		#print(pairs);
		#print(nops(pairs));
		Rprime:=subsop(i=NULL,base(R));
		ans:=ans and evalb(nops(pairs)=nops(pos_roots(Rprime))+rank(Rprime));
		for pair in pairs do
			for t in pair do
				ans:=ans and not(member(i,t));
			od;
		od;
	od;
	return ans;
end proc;

# checks condition (ii)
check_compat_ii:=proc(R,cox) local ans,i;
	ans:=true;
	for i from 1 to rank(R) do
		if nops(reduce([i,op(cox)],R))<nops(cox) then
			print(`compat (ii)`,i);
			#print([op(compat_rel_cl(R,cox))]);
			#print(map(x->map(sigma,x,R,i),[op(compat_rel_cl(R,reduce([i,op(cox),i],R)))]));
			ans:=ans and evalb(compat_rel_cl(R,cox)=map(x->map(sigma,x,R,i),compat_rel_cl(R,reduce([i,op(cox),i],R))));
		end if;
	od;
	return ans;
end proc;

# checks Proposition 7.3 {restrict}
# for each choice of coxeter element and returns
# the list of all different results.
check_compat_res_all:=proc(R) local f, num_checked,results,inttemp;
	results:=true;
	num_checked:=0;
	inttemp:=interface(quiet);
	interface(quiet=true);
	f:=proc(cox)  local restmp,t;
		t:=time();
		restmp:=check_compat_res(R,cox);
		results:=results and restmp; 
		num_checked:=num_checked+1;
		printf("%a    %a    %a    %a       %a\n",name_of(R),num_checked,cox,restmp,time()-t);
	end proc;
	f_on_cox(R,f);
	interface(quiet=inttemp);
	return(results);
end proc;

#  This checks Proposition {restrict}
#  It is sufficient to check the case where J has one fewer
#  element than S.
check_compat_res:=proc(R,cox) local ans,i,rel,relprime;
	ans:=true;
	for i from 1 to rank(R) do
		rel:=compat_rel_cl(R,cox);
		#print(rel);
		relprime:=compat_rel_cl(subsop(i=NULL,base(R)),res_cox(cox,i));
		#print(relprime);
		ans:=ans and evalb(restrict_pairs(rel,i)=shift_indices_pm(rank(R),i,relprime));
		#print(restrict_pairs(rel,i));
		#print(shift_indices_pm(rank(R),i,relprime));
	od;
	return(ans);
end proc;



#  a procedure to restrict the compatibility relation to
#  pairs which lie in the parabolic <i>.
restrict_pairs:=proc(rel,i) local pair,res_rel,add_it;
	res_rel:={};
	for pair in rel do
		#print(pair);
		add_it:=not(member(i,op(1,pair))) and not(member(-i,op(1,pair)));
		add_it:=add_it and not(member(i,op(2,pair))) and not(member(-i,op(2,pair)));
		if add_it then
			res_rel:=res_rel union {pair};
		end if;
	od;
	return(res_rel);
end proc;

#  shifts the indices in the parabolic <i>
#  so that they agree with the indices in the big group.
#  n is the rank of the big group
#  obj is any object whose operands are indices with respect 
# to the parabolic.
shift_indices:=proc(n,i,obj) local newobj,j;
	newobj:=obj;
	for j from n-1 to i by -1 do
		newobj:=subs(j=j+1,newobj);
	od;
	return newobj;
end proc;


#  shifts the indices in the parabolic <i>
#  so that they agree with the indices in the big group.
#  n is the rank of the big group
#  obj is any object whose operands are indices with respect 
# to the parabolic, or negatives of these indices.
shift_indices_pm:=proc(n,i,obj) local newobj,j;
	newobj:=obj;
	for j from n-1 to i by -1 do
		newobj:=subs(j=j+1,-j=-j-1,newobj);
	od;
	return newobj;
end proc;





# applies sigma_i to an "almost positive reflection" t
# Here an almost positive reflection is either a reduced word
#  for a reflection or a list with one element:  the negation 
#  of an index for a simple reflection.
sigma:=proc(t,R,i) local j;
	if nops(t)= 1 then
		j:=op(1,t);
		if abs(j)=i then return([-j]);
		elif j<0 then return t;
		end if;
	end if;
	return(reduce([i,op(t),i],R));
end proc;


#  this runs the function f on every sortable element.
#  R is a root system datum, cox is a word for a Coxeter element
#  (as a list of indices of simple generators of base(R)) and 
#  f is a function that either takes as an argument a c-sorting word
#  for a sortable element, or takes two arguments, first a c-sorting
#  word, and second a vector.
#  This vector is the image of interior_pt(R) under the 
#  sortable element, so that right inversions of the sortable
#  element can be determined by taking inner products (iprod)
#  of this vector with positive roots.
#  calling sequence:
#  f_on_sortas(R,cox,f)
f_on_sortas:=proc(R,cox,f)
		if crys then 
			f_o_s_helper(base(R),cox,[],interior_pt(R),f);
		else
			f_o_s_non_crys_helper(base(R),cox,[],interior_pt(R),f);
		end if;
end proc;

f_o_s_helper:=proc(R,cox,so_far,vec,f)
	if cox=[] then f(so_far,vec);
	else
		if iprod(op(op(1,cox),R),vec)>0 then
			f_o_s_helper(R,[op(2..nops(cox),cox),op(1,cox)],[op(so_far),op(1,cox)],reflect(op(op(1,cox),R),vec),f);
		end if;
		f_o_s_helper(R,[op(2..nops(cox),cox)],so_far,vec,f);
	end if;
end proc;

f_o_s_non_crys_helper:=proc(R,cox,so_far,vec,f)
	if cox=[] then f(so_far,vec);
	else
		if fuzzy_pos(iprod(op(op(1,cox),R),vec)) then
			f_o_s_non_crys_helper(R,[op(2..nops(cox),cox),op(1,cox)],[op(so_far),op(1,cox)],reflect(op(op(1,cox),R),vec),f);
		end if;
		f_o_s_non_crys_helper(R,[op(2..nops(cox),cox)],so_far,vec,f);
	end if;
end proc;

#  this runs the function f on every sortable element,
#  but here f takes two or three arguments, first the distinguished 
#  word for the sortable element and second a list of the
#  images of the simple roots under the inverse of the element.
#  the function f is also presented with a third argument, a vector
#  as explained before f_on_sortas.
#  The significance of the base_images is that s is a cover reflection of
#  w if and only if the image, under w^{-1} of the corresponding
#  simple root is the negative of a simple root.

f_on_sortas_and_base_images:=proc(R,cox,f)
	if crys then
		fosabi_helper(base(R),cox,[],interior_pt(R),base(R),f);
	else
		fosabi_noncrys_helper(base(R),cox,[],interior_pt(R),base(R),f);
	end if;
end proc;

fosabi_helper:=proc(R,cox,so_far,vec,base_images,f) local new_base_images;
	if cox=[] then f(so_far,base_images,vec);
	else
		if iprod(op(op(1,cox),R),vec)>0 then
			new_base_images:=map2(reflect,op(op(1,cox),R),base_images);
			fosabi_helper(R,[op(2..nops(cox),cox),op(1,cox)],
				[op(so_far),op(1,cox)],reflect(op(op(1,cox),R),vec),new_base_images,f);
		end if;
		fosabi_helper(R,[op(2..nops(cox),cox)],so_far,vec,base_images,f);
	end if;
end proc;

fosabi_noncrys_helper:=proc(R,cox,so_far,vec,base_images,f) local new_base_images;
	if cox=[] then f(so_far,base_images,vec);
	else
		if fuzzy_pos(iprod(op(op(1,cox),R),vec)) then
			new_base_images:=map2(reflect,op(op(1,cox),R),base_images);
			fosabi_noncrys_helper(R,[op(2..nops(cox),cox),op(1,cox)],
				[op(so_far),op(1,cox)],reflect(op(op(1,cox),R),vec),new_base_images,f);
		end if;
		fosabi_noncrys_helper(R,[op(2..nops(cox),cox)],so_far,vec,base_images,f);
	end if;
end proc;

#  this runs the function f on every sortable element,
#  but here f can take two or three arguments.
#  the first is a word as in f_on_sortas and the 
#  second is a list of positive roots corresponding to the left inversions
#  of the sortable element.
#  The function is also presented a third argument, a vector as 
#  in f_on_sortas.
#  f_on_sortas_and_inv(R,cox,f)
f_on_sortas_and_inv:=proc(R,cox,f)
	if crys then
		fosai_helper(base(R),cox,[],interior_pt(R),[],f);
	else
		fosai_non_crys_helper(base(R),cox,[],interior_pt(R),[],f);
	end if;
end proc;

fosai_helper:=proc(R,cox,so_far,vec,inv,f)
	if cox=[] then f(so_far,inv,vec);
	else
		if iprod(op(op(1,cox),R),vec)>0 then
			fosai_helper(R,[op(2..nops(cox),cox),op(1,cox)],[op(so_far),op(1,cox)],
					reflect(op(op(1,cox),R),vec),[op(inv),reflect(seq(op(i,R),i=so_far),op(op(1,cox),R))],f);
		end if;
		fosai_helper(R,[op(2..nops(cox),cox)],so_far,vec,inv,f);
	end if;
end proc;

fosai_non_crys_helper:=proc(R,cox,so_far,vec,inv,f)
	if cox=[] then f(so_far,inv,vec);
	else
		if fuzzy_pos(iprod(op(op(1,cox),R),vec)) then
			fosai_non_crys_helper(R,[op(2..nops(cox),cox),op(1,cox)],[op(so_far),op(1,cox)],
					reflect(op(op(1,cox),R),vec),[op(inv),reflect(seq(op(i,R),i=so_far),op(op(1,cox),R))],f);
		end if;
		fosai_non_crys_helper(R,[op(2..nops(cox),cox)],so_far,vec,inv,f);
	end if;
end proc;

# f_on_cox apples a function f to all different group
# elements which can be obtained as a product of
# simple generators (each generator occurring once).
# the elements are output as words in the indices of the generators
# ie the standard format for reduced words in the coxeter package.
# WARNING: should be applied only to irreducible groups
# the reducible case could easily be programmed, but I had no need. 
# might also give wierd results for the rank zero Coxeter group...
f_on_cox:=proc(R,f) local i,j,mat,orient, antiorient,numedges,f_for_subs;
	mat:=cox_matrix(R);
	orient:=[];
	antiorient:=[];
	for i from 1 to rank(R) do
		for j from i+1 to rank(R) do
			if mat[i,j]>2 then 
				orient:=[op(orient),[i,j]];
				antiorient:=[op(antiorient),[j,i]];
			end if;
		od;
	od;
	#print(orient,antiorient);
	numedges:=nops(orient);
# this will be a function applied to each sublist of the list [1..numedges]
	f_for_subs:=proc(list) local k,poset;
		poset:={};
		for k from 1 to numedges do
			if k in list then 
				poset:=poset union {op(k,orient)};
			else
				poset:=poset union {op(k,antiorient)};
			end if;
		od;
		#print(list,poset);
		#print(extensions(poset));
		f(op(1,extensions(poset)));
	end proc;
	f_on_sublists([seq(i,i=1..numedges)],f_for_subs);
end proc;


#  This finds the cox-sorting word for the longest element w_0.
#  R is a root system datum, cox is a word for a Coxeter element
#  (as a list of indices of simple generators of base(R)) and
#  output will be a cox-sorting word for the longest element 
#  (again as a list of indices of simple generators of base(R))
#  calling sequence:
#  w0sorting(R,cox)
w0sorting:=proc(R,cox)
	if crys then
		w0sort_helper(base(R),cox,[],interior_pt(R),interior_pt(R));
	else
		w0sort_non_crys_helper(base(R),cox,[],interior_pt(R),interior_pt(R));
	end if;
end proc;

w0sort_helper:=proc(R,cox,so_far,vec,int)
	if -vec=int then
		return(so_far);
	elif iprod(op(op(1,cox),R),vec)>0 then
		return(w0sort_helper(R,[op(2..nops(cox),cox),op(1,cox)],[op(so_far),op(1,cox)],reflect(op(op(1,cox),R),vec),int));
	else	
		return(w0sort_helper(R,[op(2..nops(cox),cox),op(1,cox)],so_far,vec,int));
	end if;
end proc;

w0sort_non_crys_helper:=proc(R,cox,so_far,vec,int)
	if fuzzy_vec_eq(-vec,int) then
		return(so_far);
	elif fuzzy_pos(iprod(op(op(1,cox),R),vec)) then
		return(w0sort_non_crys_helper(R,[op(2..nops(cox),cox),op(1,cox)],[op(so_far),op(1,cox)],reflect(op(op(1,cox),R),vec),int));
	else	
		return(w0sort_non_crys_helper(R,[op(2..nops(cox),cox),op(1,cox)],so_far,vec,int));
	end if;
end proc;

# gives the reflection ordering on positive roots associated to the 
# cox sorting word for the longest element.
# calling sequence:
# ref_order(R,cox) 
ref_order:=proc(R,cox) local order,word,i; 
	order:=[];
	word:=w0sorting(R,cox);
	for i from 1 to nops(word) do
		order:=[op(order),reflect(seq(op(op(j,word),base(R)),j=1..i))];		
	od;
	return order;
end proc;






#  A helper function which restricts a coxeter element
#  (a list of indices in the base) to a coxeter element
#  of the parabolic obtained by deleting the ith entry of the base.
res_cox:=proc(cox,i) local new_cox,j;
	new_cox:=subs(i=NULL,cox);
	#print(new_cox);
	for j from 1 to nops(new_cox) do
		if op(j,new_cox)> i then 
			new_cox:=subsop(j=op(j,new_cox)-1,new_cox);
			#print(new_cox);
		end if;
	od;
	return(new_cox);
end proc;

#  This applies the map cl to any sortable element.
#  Output is a set of "almost positive reflections."
#  The input is assumed to be the c-sorting word for 
#  the element, so the procedure  does not explicitly need 
#  to know what c is.
cl:=proc(R,sorta) local last_refs,i,j,l,in_supp;
	last_refs:={};
	in_supp:=array(1..rank(R));
	for i from 1 to rank(R) do
		in_supp[i]:=false;
	od;
	for j from 1 to nops(sorta) do
		in_supp[op(j,sorta)]:=true;
		if not(member(op(j,sorta),[op(j+1..nops(sorta),sorta)])) then
			last_refs:=last_refs union {reduce([op(1..j,sorta),seq(op(j-k,sorta),k=1..j-1)],R)};
		end if;
	od;
	for l from 1 to rank(R) do
		if not(in_supp[l]) then
			last_refs:=last_refs union {[-l]};
		end if;
	od;
	return(last_refs);
end proc;

#  This applies the map cl to any sortable element.
#  Output is a set of roots.
#  The input is assumed to be the c-sorting word for 
#  the element, so the procedure  does not explicitly need 
#  to know what c is.
cl_roots:=proc(R,sorta) local B,last_roots,i,j,l,in_supp;
	B:=base(R);
	last_roots:={};
	in_supp:=array(1..rank(R));
	for i from 1 to rank(R) do
		in_supp[i]:=false;
	od;
	for j from 1 to nops(sorta) do
		in_supp[op(j,sorta)]:=true;
		if not(member(op(j,sorta),[op(j+1..nops(sorta),sorta)])) then
			last_roots:=last_roots union {reflect(seq(op(op(m,sorta),B),m=1..j))};
		end if;
	od;
	for l from 1 to rank(R) do
		if not(in_supp[l]) then
			last_roots:=last_roots union {-op(l,B)};
		end if;
	od;
	return(last_roots);
end proc;

#  This computes the compatibility relation based on the 
#  following definition:  two distinct "almost positive reflections" 
#  are compatible if and only if there is some c-sortable element 
#  w such that both are in cl_c(w).
#  the output is a set of 2-element sets, the collection of all compatible pairs.
compat_rel_cl:=proc(R,cox) local rel,f;  option remember;
	rel:={};
	f:=proc(sorta) local last_refs,i,j;
		last_refs:=[op(cl(R,sorta))];
		for i from 1 to nops(last_refs) do
			for j from i+1 to nops(last_refs) do
				rel:=rel union{{op(i,last_refs),op(j,last_refs)}};
			od;
		od;
	end proc;
	f_on_sortas(R,cox,f);
	return rel;
end proc;

# orthogonality of vectors (exact)
orth:=proc(a,b)
	return(evalb(iprod(a,b)=0));
end;

#  subset test, up to the coxeter package's epsilon
#  elements of the sets are vectors
fuzzy_vec_subset:=proc(set1,set2) local list1,list2,x,y;
	list1:=[op(set1)];
	list2:=[op(set2)];
	for x in list1 do
		if not(fuzzy_vec_member(x,list2)) then
			return(false);
		end if;
	od;
	return true;
end proc;

#  Tests set equality of sets of vectors, up to epsilon
#  NOTE:  This is not the fastest possible way to check this,
#  but it's good enough for my purposes.
fuzzy_vec_set_eq:=proc(set1,set2)
	return(fuzzy_vec_subset(set1,set2) and fuzzy_vec_subset(set2,set2));
end proc;



# membership in a list up to the coxeter package's epsilon.
# just true/false.  
# also has the position option as in the member command:
#  if a third argument is given, it will be assigned to the first position
# of the vector in the list.
fuzzy_vec_member:=proc(a,list) local i;
	for i from 1 to nops(list) do
		if fuzzy_vec_eq(op(i,list),a) then 
			if nargs=3 then assign(args[3],i); end if;
			return true; 
		end if;
	od;
	return false;
end;

default_squared:=`coxeter/default`[epsilon]^2;

# equality of vectors up to the coxeter package's epsilon.
fuzzy_vec_eq:=proc(a,b)
	return(evalb(iprod(a-b,a-b)<default_squared));
end;

# equality of scalars to the coxeter package's epsilon.
fuzzy_eq:=proc(a,b)
	return(fuzzy_zero(a-b));
end;

# orthogonality of vectors up to the epsilon
fuzzy_orth:=proc(a,b)
	return(fuzzy_zero(iprod(a,b)));	
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;

# 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;

#  checks whether the two given elements or their opposites occur in the given 
#  list in the given order;
pm_in_order:=proc(x,y,list) local len,found_x,i,j,goodi;
	len:=nops(list);
	found_x:=false;
	for i from 1 to len do
		if op(i,list)=x or op(i,list)=-x then
			goodi:=i;
			found_x:=true;
			break;
	      	end if;
	od;
	if found_x then
		for j from goodi+1 to len do
		   	if op(j,list)=y or op(j,list)=-y then
				return(true);
			end if;
		od;
	end if;
	return(false);
end proc;

#  a fuzzy version.  For Vectors only!
fuzzy_pm_in_order:=proc(x,y,list) local len,found_x,i,j,goodi;
	len:=nops(list);
	found_x:=false;
	for i from 1 to len do
		if fuzzy_vec_eq(op(i,list),x) or fuzzy_vec_eq(op(i,list),-x) then
			goodi:=i;
			found_x:=true;
			break;
	      	end if;
	od;
	if found_x then
		for j from goodi+1 to len do
		   	if fuzzy_vec_eq(op(j,list),y) or fuzzy_vec_eq(op(j,list),-y) then
				return(true);
			end if;
		od;
	end if;
	return(false);
end proc;


#  returns false if x and y (vectors) are both non-orthogonal and
#  out of order in list.
#  Returns true otherwise.
#  assumes there are no repeats in list.
not_misordered:=proc(x,y,list) local i,j;
	print(iprod(x,y));
	if not(orth(x,y)) then
		for i from 1 to nops(list) do
			if op(i,list)=y then
				print(i);
				for j from i+1 to nops(list) do
					if op(j,list)=x then
						return(false);
					end if;
				od;
				break;
			end if;
		od;
	end if;
	return(true);
end proc;

#  a fuzzy version
fuzzy_not_misordered:=proc(x,y,list) local i,j;
	print(iprod(x,y));
	if not(fuzzy_orth(x,y)) then
		for i from 1 to nops(list) do
			if fuzzy_vec_eq(op(i,list),y) then
				print(i);
				for j from i+1 to nops(list) do
					if fuzzy_vec_eq(op(j,list),x) then
						return(false);
					end if;
				od;
				break;
			end if;
		od;
	end if;
	return(true);
end proc;

# finds those indices which shorten cox on the left
firsts:=proc(R,cox) local i,f;
	f:=[];
	for i from 1 to rank(R) do
		if nops(reduce([i,op(cox)],R))<nops(cox) then
			f:=[op(f),i];
		end if;
	od;
	return(f);
end proc;

# finds those indices which shorten cox on the right
lasts:=proc(R,cox) local i,l;
	l:=[];
	for i from 1 to rank(R) do
		if nops(reduce([op(cox),i],R))<nops(cox) then
			l:=[op(l),i];
		end if;
	od;
	return(l);
end proc;

#  tests if a root is positive
is_positive:=proc(R,r) local co,i;
	co:=root_coords(r,R);
	if crys then
		for i from 1 to rank(R) do
			if op(i,co)>0 then return(true);
			elif op(i,co)<0 then return(false);
			end if;
		od;
	else
		for i from 1 to rank(R) do
			if fuzzy_pos(op(i,co)) then return(true);
			elif fuzzy_neg(op(i,co)) then return(false);
			end if;

		od;

	end if;
end proc;

# Given a reduced word for w and 
# a vector representing the image of 
# interior_pt(R) (or some other point in the interior
# of the fundamental chamber) under w,
# this finds the set of cover reflections of w,
# as positive roots.
# R is the root system datum
cov_refs:=proc(R,word,vec) local b, cov, i;
	b:=base(R);
	cov:={};
	if crys then
		for i from 1 to nops(b) do
			if iprod(op(i,b),vec)<0 then
				cov:=cov union {reflect(seq(op(j,b),j=word),op(i,b))};
			end if;
		od;
	else
		for i from 1 to nops(b) do
			if fuzzy_neg(iprod(op(i,b),vec)) then
				cov:=cov union {reflect(seq(op(j,b),j=word),op(i,b))};
			end if;
		od;
	end if;
	return(cov);
end proc;

#  tests the cov_refs procedure by calling it on every c-sortable element
test_cov_refs:=proc(R,cox) local f;
	f:=proc(word,vec)
		print(word,cov_refs(R,word,vec));
	end proc;
	f_on_sortas(R,cox,f);
end proc;

#  gets R,cox and pr, where pr is pos_roots(R).
#  Returns a pair ji_inv_sets,mi_inv_sets
#  ji_inv_sets is an array with entries counted by the positive roots.
#  each entry i is inversion set of the join-irreducible c-sortable 
#  element whose unique cover reflection is the ith entry of pr.
#  mi_inv_sets is an array whose ith entry is the inversion set
#  of the meet-irreducible c-sortable element forming an 
#  interval [v,w] (see discussion of Conjecture {fan iso})
#  with the ith entry of ji_inv_sets.
#  Exception:  if the w corresponding to a join-irreducible
#  is w_0, then the corresponding entry of mi_inv_sets
#  is the symbol 'w0'  
find_intervals:=proc(R,cox,pr) local B,num_refs,jis,mis,pds,mi_inv_sets_misordered,f,ji_inv_sets,mi_inv_sets,k,l;
	B:=base(R);
	num_refs:=nops(pr);
	jis:=array[1..numrefs];
	ji_inv_sets:=array[1..numrefs];
	mis:=[];
	pds:=[];
	mi_inv_sets_misordered:=[];
	if crys then
		f:=proc(word,inv,vec) local i,j,rt_des,last_vec,k,pd;
			rt_des:=[];
			for i from 1 to rank(R) do
				if iprod(op(i,B),vec)<0 then
					rt_des:=[op(rt_des),i];
				end if;
			od;
			#print(word,rt_des);
			if nops(rt_des)=1 then 
				last_vec:=op(nops(inv),inv);
				member(last_vec,pr,'j');
				jis[j]:=word;
				ji_inv_sets[j]:=inv;
			end if;
			if nops(rt_des)=rank(R)-1 then
				k:=op(1,{seq(z,z=1..rank(R))} minus {op(rt_des)});
				#print(k);
				mis:=[word,op(mis)];
				mi_inv_sets_misordered:=[inv,op(mi_inv_sets_misordered)];
				pd:=pidown(R,cox,
					reduce([op(word),op(shift_indices(rank(R),k,longest_elt([seq(op(m,B),m=rt_des)])))],R)
				);
				#print(word,rt_des,shift_indices(rank(R),k,longest_elt([seq(op(m,B),m=rt_des)])),pd);
				pds:=[pd,op(pds)];
			end if;
		end proc;
	else
		f:=proc(word,inv,vec) local i,j,rt_des,last_vec,k,pd;
			rt_des:=[];
			for i from 1 to rank(R) do
				if fuzzy_neg(iprod(op(i,B),vec)) then
					rt_des:=[op(rt_des),i];
				end if;
			od;
			#print(word,rt_des);
			if nops(rt_des)=1 then 
				last_vec:=op(nops(inv),inv);
				fuzzy_vec_member(last_vec,pr,'j');
				jis[j]:=word;
				ji_inv_sets[j]:=inv;
			end if;
			if nops(rt_des)=rank(R)-1 then
				k:=op(1,{seq(z,z=1..rank(R))} minus {op(rt_des)});
				#print(k);
				mis:=[word,op(mis)];
				mi_inv_sets_misordered:=[inv,op(mi_inv_sets_misordered)];
				pd:=pidown(R,cox,
					reduce([op(word),op(shift_indices(rank(R),k,longest_elt([seq(op(m,B),m=rt_des)])))],R)
				);
				#print(word,rt_des,longest_elt([seq(op(m,B),m=rt_des)]),pd);
				pds:=[pd,op(pds)];
			end if;
		end proc;
	end if;
	f_on_sortas_and_inv(R,cox,f);
	#print(jis,ji_inv_sets,mis,mi_inv_sets_misordered,pds);
	mi_inv_sets:=array[1..num_refs];
	for k from 1 to num_refs do mi_inv_sets[k]:='w0' od;
	for k from 1 to num_refs do
		for l from 1 to num_refs do
			#print(jis[l],op(k,pds),op(k,mis),evalb(jis[l]=op(k,pds)));
			if jis[l]=op(k,pds) then
				mi_inv_sets[l]:=op(k,mi_inv_sets_misordered);
				break;
			end if;
		od;
	od;	
	#  for the remaining join-irreducibles not hit by the above for-loop,
	#  the symbol 'w0' remains.
	#  this is valid because each of these join-irreducibles
	#  must be $\pidown(w(w_0)_J)$ for some pair $(w,J)$ 
	#  if some join-irreducible wasn't hit by a meet-irreducible
	# in the above fo loop, the pair $(w,J)$ must have $w=w_0$
	#print(ji_inv_sets,mi_inv_sets);
	return(ji_inv_sets,mi_inv_sets);
end proc;

#  does the map pidown
#  starting with a reduced word, gives the cox-sorting word.
#  If the input word is not reduced, this could give a wierd answer.
#  This is recursive, and in lower levels of the recursion, cox is 
#  a coxeter element for a standard parabolic subgroup.
pidown:=proc(R,cox,word) local sw;
	if cox=[] then return []; end if;
	#print([op(1,cox),word],R);
	sw:=reduce([op(1,cox),op(word)],R);
	if nops(sw)<nops(word) then
		return([op(1,cox),op(pidown(R,[seq(op(i,cox),i=2..nops(cox)),op(1,cox)],sw))]);
	else
		return(pidown(R,[seq(op(i,cox),i=2..nops(cox))],word));
	end if;
end proc;