# ---------------------------------------------------------------
# Code to compute Schubert calculus using the recurrence from
# Allen Knutson's paper "A Schubert calculus recurrence from
# the noncomplex W-action on G/B".
#
# This code was written, in part, for an upcoming joint project
# with Hugh Thomas, I greatly acknowledge Hugh for motivating
# conversations. 
# 
# Author: Alexander Yong
# report errors to Alexander Yong ayong@math.umn.edu
# Date: June 4, 2007
# v0.2: changes the code slightly from Knutson's paper, to correct for a coefficient
# error; basically, (\alpha,\beta) should be (\alpha^{\vee},\beta).
#
# Uses: Stembridge's coxeter2_2v.txt, traverse, software
#
# Quick instructions:
#
# 1. Get coxeter2_2v.txt and traverse.txt from John
#    Stembridge's website
#    http://www.math.lsa.umich.edu/~jrs
# 2. Start maple
# 3. type in ``read "../coxeter2_2v.txt";'' and ``withcoxeter();''
# 4. type in ``read "../traverse.txt";''
# 5. type in ``read "../Schubert.v0.1.txt";''
# 6. ``Lie type'' refers to A1, A2, A3, ..., B2,...,E6,..., see
#    coxeter2_2v.txt's conventions
# 7. Elements of the Weyl group are interpreted as
#    words in the simple relfections, e.g., [1,2,1]=[2,1,2] in
#    Lie type A3. Again, see coxeter2_2v.txt's conventions
# 8. "compute_single_SSC([1],[2], [2,1], E6)" computes C_{[1], [2]}^{[2,1]}
#    in type E6
# 9. "Schub_product([1],[2], B5)" compute the product of Schubert
#10. classes X[1]*X[2] in the basis of Schubert classes in type B5.
#11. The Schub_product algorithm is pretty slow. For type A, there is much
#    faster software available (ACE for Maple).  
# ---------------------------------------------------------------


# Recursive code to compute a single Schubert structure constant
# c_{wv}^{u}

`compute_single_SSC`:=proc(w,v,u, Lie_type) 
	local ii, answer, flag, simple_found_flag, w_length, alpha, big_term,
              jj, w_prime, inn_coeff;

	flag:=0:
	if(reduce([op(w), op(longest_elt(Lie_type))], Lie_type)=[]) then
		flag:=1:
		if(reduce([op(u), op(longest_elt(Lie_type))], Lie_type)=[]) then
			if(reduce(v, Lie_type)=[]) then
				answer:=1:
			else
				answer:=0:
			fi:
		else
			answer:=0:
		fi:
		return(answer):
	fi:

	if(flag=0) then
		# find a simple reflection that increases the length (by 1)
		ii:=1:
		simple_found_flag:=0:
		w_length:=nops(reduce(w, Lie_type)):
		while(simple_found_flag=0 and ii<=rank(Lie_type)) do
			if(nops(reduce([op(w), ii], Lie_type))>w_length) then
				simple_found_flag:=1:
				alpha:=ii:
			else
				ii:=ii+1:
			fi:
		od:

		# check dc-triviality
		if ((nops(reduce([op(v), alpha], Lie_type))>nops(reduce(v, Lie_type))) and
		    (nops(reduce([op(u), alpha], Lie_type))<nops(reduce(u, Lie_type)))) then
		    answer:=0:
		    return(answer):
		    flag:=1:
		fi:

		if(flag=0) then
		    if ((nops(reduce([op(v), alpha], Lie_type))>nops(reduce(v, Lie_type))) and
		    (nops(reduce([op(u), alpha], Lie_type))<nops(reduce(u, Lie_type)))) then 

			answer:=compute_single_SSC(reduce([op(w), alpha], Lie_type), reduce([op(v),alpha], Lie_type), 
                                reduce(u, Lie_type), Lie_type):
			return(answer):
		    elif ((nops(reduce([op(v), alpha], Lie_type))>nops(reduce(v, Lie_type))) and
		    (nops(reduce([op(u), alpha], Lie_type))>nops(reduce(u, Lie_type)))) then
			
                        answer:=compute_single_SSC(reduce([op(w), alpha], Lie_type), reduce(v, Lie_type), 
                                reduce([op(u),alpha], Lie_type), Lie_type):
			RETURN(answer): 		    
		    else
		        
		        big_term:=0:
			w_length:=nops(reduce(w, Lie_type)):
			for jj from 1 to nops(pos_roots(Lie_type)) do
			    if(op(jj,pos_roots(Lie_type))<>op(alpha, base(Lie_type))) then
					
					w_prime:=reduce([op(w), op(pos_root_to_word(jj, Lie_type))], Lie_type):
					
					if(nops(reduce(w_prime, Lie_type))=(w_length+1)) then
								
						inn_coeff:=iprod((coify(op(jj,pos_roots(Lie_type))),(op(alpha,base(Lie_type))))):					
						#inn_coeff:=simplify((op(jj,pos_roots(Lie_type))
						#-reflect((op(alpha,base(Lie_type))),
			                        #           (op(jj,pos_roots(Lie_type)))))/(op(alpha,base(Lie_type)))):
						big_term:=big_term+inn_coeff*compute_single_SSC(w_prime, 
                                                          reduce([op(v), alpha], Lie_type), reduce(u, Lie_type), Lie_type):
					fi:
			    fi:
			od:
			
			answer:=compute_single_SSC([op(w), alpha], v, [op(u),alpha], Lie_type)
                                +compute_single_SSC(reduce([op(w), alpha], Lie_type), reduce([op(v),alpha], Lie_type), 
                                 reduce(u, Lie_type), Lie_type)+big_term:			
			return(answer):
		    fi:
		fi:
	fi:
end:

`coify`:=proc(avector)
	local newavector;
	#print(avector);
	if(avector<>0) then
		       newavector:=2*avector/(iprod(avector,avector)):
	else
		       newavector:=0:
        fi:		       
	RETURN(newavector):
end:

`pos_root_to_word`:=proc(index_pos_root, Lie_type)
	 local c0, r, w, v0;
	 v0:=interior_pt(Lie_type):
	 r:=op(index_pos_root, pos_roots(Lie_type)):
	 vec2fc(reflect(r,v0), Lie_type, 'w'): 
	 RETURN(w):
end:

# Code to multiply two Schubert classes in the generalized
# flag manifold of a given Lie type
# This procedure is probably very slow, and is used
# when there is no extra information known about the
# product taken.
# Modifications of this procedure will include improvements
# in performance coming from that inhanced data

`Schub_product`:=proc(w,v, Lie_type)
	local ii, w_length,  v_length, u_length, the_list, the_poly;
	
	# first build the list of all permutations of the
	# appropriate length

	w_length:=nops(reduce(w, Lie_type)):
	v_length:=nops(reduce(v, Lie_type)):	
	u_length:=w_length+v_length:
	the_list:=traverse(Lie_type,[],length=u_length,find_length_elts, u_length):
	the_poly:=0:
	for ii from 1 to nops(the_list) do
	    the_poly:=the_poly+compute_single_SSC(w,v, op(ii, the_list), Lie_type)*X[op(ii,the_list)]:	    
	    
	od:
	return(the_poly):
end:

# use traverse to find elements of a given length
# call by "traverse(Lie_type, [], find_length_elts, length):"

`find_length_elts`:=proc(list_so_far, w,v, Lie_type, length_wanted)
	if(nops(reduce(w, Lie_type))=length_wanted) then
		return([op(list_so_far), w]):
	else
		return(list_so_far):
	fi:
end: