with(coxeter):


#  To calculate f-vector correction factors for the exceptional
#  types, do:

# read(f_maple):
#for R in [E6,E7,E8,F4,H3,H4] do
#	printf("\n");
#	print(R,cox_number(R));
#	for i from 0 to rank(R) do 
#			print(i,correction_f(R,i)); 
#	od;
#od;

# For the h-vector correction factors, replace correction_f
# by correction_h

# takes a root system data structure (of finite type) and returns
# a polynomial such that the coeff of x^(n-k) is f_k.
f_poly:=proc(R) option remember;
	return(x^rank(R)*subs(x=1/x,f_poly_rev(R)));
end;

# takes a root system data structure (of finite type) and returns
# reversed from what is defined above.
# that is, the coeff of x^k is f_k.
f_poly_rev:=proc(R) local name,first; option remember;
	if rank(R)=0 then return 1; end if;
	if irreduc(name_of(R)) then return f_poly_irr_rev(R);
	else
		name:=name_of(R);
		first:=op(1,op(2,factors(name)));
		return (f_poly_rev(op(1,first))^op(2,first))*f_poly_rev(name/op(1,first)^op(2,first));
	end if;
end;

# works on a finite irreducible Coxeter group.
f_poly_irr_rev:=proc(R) local partial_f_poly,f; option remember;
	partial_f_poly:=0;
	f:=proc(x) ;
		partial_f_poly:=partial_f_poly+f_poly_rev(x);
	end proc; 
	f_on_k_sublists(base(R),rank(R)-1,f);
	return((cox_number(R)*m+2)*subs(xprime=x,int(partial_f_poly,x=0..xprime,'continuous'))/2+1);
end;

# the reciprocal f_polynomial
recip_f_poly:=proc(R)
	return(subs(x=-x,m=-m-1,f_poly(R)));
end;

# an entry in the h-vector
h:=proc(R,k)
	return(factor(coeff(factor(h_poly(R)),x,rank(R)-k)));
end;

# an entry in the f-vector
# recall that f_k is the number of k-element
#simplices, not k-dimensional simplices!
f:=proc(R,k)
	return(factor(coeff(factor(f_poly(R)),x,rank(R)-k)));
end;

# takes a root system data structure (of finite type) and returns
# what we hope will be the associated Narayana polynomial.
h_poly:=proc(R) option remember;
	return subs(x=x-1,f_poly(R));	
end;

# takes a root system data structure (of finite type) and returns
# what we hope will be the associated Narayana polynomial.
recip_h_poly:=proc(R) option remember;
	return subs(x=x-1,recip_f_poly(R));	
end;

# an entry in the reciprocal h-vector
recip_h:=proc(R,k)
	return(factor(coeff(factor(recip_h_poly(R)),x,rank(R)-k)));
end;

# lists the exponents of level =< k.
level_le_k:=proc(R,k) local letter,n;
	if k=0 then return([]); end if;
	letter:=substring(R,1..1);
	n:=convert(substring(R,2..length(R)),decimal,10);
	if letter=D then
		if k=1 or k=2 then return([1]);
		elif k<n-1 then return([seq(2*i-1,i=1..k-1)]);
		elif k=n-1 then return(sort([n-1,seq(2*i-1,i=1..n-2)]));
		elif k=n then return(sort([n-1,seq(2*i-1,i=1..n-1)]));
		end if;
	elif R in [E6,E7,F4,H4] then
		if k=1 or k=2 then return([op(1,exponents(R))]);
		elif k<n-1 then return([op(1..k-1,exponents(R))]);
		elif k=n-1 then return([op(1..k,exponents(R))]);
		elif k=n then return(exponents(R));
		end if;
	elif R=E8 then
		if k=1 or k=2 then return([op(1,exponents(R))]);
		elif k=3 or k=4 then return([op(1..2,exponents(R))]);
		elif k=5 or k=6 then return([op(1..k-1,exponents(R))]);
		elif k=7 then return([op(1..7,exponents(R))]);
		elif k=8 then return(exponents(R));
		end if;
	elif k<= n then return([op(1..k,exponents(R))]);	
	end if;
end;

# only for irreducibles!!!
# gives the correction factor to f_k
correction_f:=proc(R,k) local corr,expo; option remember;
	corr:=f(R,k);
	for expo in level_le_k(R,k) do
		corr:=corr/(m*cox_number(R)+expo+1)*(expo+1);
	od;
	corr:=corr/binomial(rank(R),k);
	return(factor(corr));	
end;

# only for irreducibles!!!
# gives the correction factor to h_k
correction_h:=proc(R,k) local corr,expo; option remember;
	corr:=h(R,k);
	for expo in level_le_k(R,k) do
		corr:=corr/(m*cox_number(R)-expo+1)*(expo+1);
	od;
	corr:=corr/binomial(rank(R),k);
	return(factor(corr));	
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;