# ---------------------------------------------------------------
# COMINRULE.v1.0.txt
# ---------------------------------------------------------------
#
# Code to compute the (co)minuscule Schubert calculus rule from
# "A combinatorial rule for cominuscule Schubert calculus", 
# coauthored by Hugh Thomas and Alexander Yong
# 
# Author: Hugh Thomas and Alexander Yong
# report errors to Alexander Yong ayong@math.umn.edu
# Date: August 6, 2006
#
# Quick instructions:
# > read ".../cominrule.v0.1.txt";
# The indexing of the code corresponds to that in the above
# mentioned paper
# cominrule(lambda, mu, nu, flagvarietytype,Lie_index, printoutput)
# is the format where
#
# flagvarietytype <-> cominuscule G/P
# -----------------------------
#        A        Grassmannian
#       LG        Lagrangian Grassmannian
#       OG        Orthogonal Grassmannian
#     oddQ        Odd quadrics
#    evenQ        Even quadrics
#        E        one of the two E-type spaces
#
# Lie_index is an integer in all cases except A, where it is
# a list [k,n].
#
# printoutput=1 causes the code to print out the SYT that
# rectify to the consecutive filling tableau (and thus
# contribute to the coefficient; printoutput=0 supresses that
#
# Examples (from the paper):
#
# cominrule([3,1],[2,1],[4,2,1],A,[4,7],1); # Gr(4,7)
# cominrule([2,1],[2,1],[4,2],LG,4,1); # LG(4,8)
# cominrule([1,1],[1,1],[1,1,1,1],oddQ,4,1); # {\mathbb Q}^{7}
# cominrule([1,1,1,2],[1,1,1,1,1],[1,1,1,2,2,1,1,1],evenQ,6,1); # {\mathbb Q}^{10}
# cominrule([1,1,1,2,5,3],[1,1,1,2,1],[1,1,1,2,5,5,2,1,1],E,7,1); #E7 
#
# and an E6 example, with output:
#
# cominrule([1,1,2,1,1],[1,1,2,2,1],[1,1,2,4,4,1],E,6,1);
#  [Z    Z    Z    4    6    _    _    _]
#  [                                    ]
#  [Z    Z    Z    2    5    7    Z    Z]
#  [                                    ],
#  [Z    Z    _    1    3    Z    Z    Z]
#  [                                    ]
#  [_    _    _    _    _    Z    Z    Z]
#
#        [Z    Z    Z    _    _    _    _    _]
#        [                                    ]
#        [Z    Z    Z    _    _    _    Z    Z]
#        [                                    ]
#        [Z    Z    4    6    _    Z    Z    Z]
#        [                                    ]
#        [1    2    3    5    7    Z    Z    Z]
#
#  [Z    Z    Z    6    7    _    _    _]
#  [                                    ]
#  [Z    Z    Z    2    4    5    Z    Z]
#  [                                    ],
#  [Z    Z    _    1    3    Z    Z    Z]
#  [                                    ]
#  [_    _    _    _    _    Z    Z    Z]
#
#        [Z    Z    Z    _    _    _    _    _]
#        [                                    ]
#        [Z    Z    Z    _    _    _    Z    Z]
#        [                                    ]
#        [Z    Z    4    6    _    Z    Z    Z]
#        [                                    ]
#        [1    2    3    5    7    Z    Z    Z]
#
#                      2
#
# Notes: We decided to do this in Maple in order for compatibility
#        of usewith John Stembridge's SF/Coxeter packages as well as 
#        the implementation of Knutson's G/B algorithm found
#        at http://www.math.umn.edu/~ayong
#        (These are not needed in this version, however.)
#
#	 This version correct an bug in v0.1 of this program,
#        where the jdt was executed incorrectly in some cases.
# ---------------------------------------------------------------

# Code to use jeu de taquin to compute a cominuscule Schubert
# structure constant. The conventions are as in the aformenentioned
# paper.

`cominrule`:=proc(lambdain,muin, nuin, flagvarietytype, Lie_index, printoutput)
	local ii, aa, bb, ambientmatrix, legalpositions, maxnumberinserted,
              lambda, mu, nu, jj, legal_positions, diff,
              lower_boundary, ans;
	
	lambda:=lambdain:
	mu:=muin:
	nu:=nuin:

	maxnumberinserted:=0;
	for ii from 1 to nops(nu) do
	    maxnumberinserted:=maxnumberinserted+op(ii,nu):	
	od:
	for ii from 1 to nops(lambda) do
	    maxnumberinserted:=maxnumberinserted-op(ii,lambda):	
	od:

	# create an array and initiate legal positions
	# which is a list of lists: 
        # ([lowest valid row positions for column 1, highest], ditto for column 2, ...)
	# based on the shape of \Lambda_{G/P} and lambda,nu
	
	if(flagvarietytype=A) then
		# Here I assume "Lie_index"=[k,n]
		ambientmatrix:=linalg[matrix](op(1,Lie_index),
                                       op(2,Lie_index)-op(1,Lie_index)):
		for ii from 1 to linalg[rowdim](ambientmatrix) do
		    for jj from 1 to linalg[coldim](ambientmatrix) do
			ambientmatrix[ii,jj]:=_:
		    od:
		od:

		# create the legal_positions matrix
		legal_positions:=linalg[matrix](2,linalg[coldim](ambientmatrix)):
		for aa from nops(lambda)+1 to linalg[coldim](ambientmatrix) do
		    lambda:=[op(lambda),0]:
		od:
		for aa from nops(nu)+1 to linalg[coldim](ambientmatrix) do
		    nu:=[op(nu),0]:
		od:
		
		diff:=nu-lambda:

		for aa from 1 to linalg[coldim](ambientmatrix) do
		    if(op(aa,diff)=0) then
			legal_positions[1,aa]:=infinity:
			legal_positions[2,aa]:=infinity:

		    else
			legal_positions[1,aa]:=op(aa,nu):
			legal_positions[2,aa]:=legal_positions[1,aa]-op(aa,diff)+1:

			legal_positions[1,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[1,aa]+1:
			legal_positions[2,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[2,aa]+1:
		    fi:
		od:
	fi:


	if(flagvarietytype=LG) then
		ambientmatrix:=linalg[matrix](Lie_index,Lie_index):
		for ii from 1 to linalg[rowdim](ambientmatrix) do
		    for jj from 1 to linalg[coldim](ambientmatrix) do
			if(jj>linalg[rowdim](ambientmatrix)-ii+1) then
				ambientmatrix[ii,jj]:=Z:
			else
				ambientmatrix[ii,jj]:=_:
			fi:
		    od:
		od:

		# create the legal_positions matrix
		legal_positions:=linalg[matrix](2,linalg[coldim](ambientmatrix)):
		for aa from nops(lambda)+1 to linalg[coldim](ambientmatrix) do
		    lambda:=[op(lambda),0]:
		od:
		for aa from nops(nu)+1 to linalg[coldim](ambientmatrix) do
		    nu:=[op(nu),0]:
		od:

		lower_boundary:=[]:
		for aa from 1 to linalg[coldim](ambientmatrix) do
		    lower_boundary:=[op(lower_boundary),aa-1]:
		od:
		
		diff:=nu-lambda:

		for aa from 1 to linalg[coldim](ambientmatrix) do
		    if(op(aa,diff)=0) then
			legal_positions[1,aa]:=infinity:
			legal_positions[2,aa]:=infinity:

		    else
			legal_positions[1,aa]:=op(aa,lower_boundary)+op(aa,nu):
			legal_positions[2,aa]:=legal_positions[1,aa]-op(aa,diff)+1:

			legal_positions[1,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[1,aa]+1:
			legal_positions[2,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[2,aa]+1:
		    fi:
		od:
	fi:

	if(flagvarietytype=OG) then
		ambientmatrix:=linalg[matrix](Lie_index-2,Lie_index-2):
		for ii from 1 to linalg[rowdim](ambientmatrix) do
		    for jj from 1 to linalg[coldim](ambientmatrix) do
			if(jj>linalg[rowdim](ambientmatrix)-ii+1) then
				ambientmatrix[ii,jj]:=Z:
			else
				ambientmatrix[ii,jj]:=_:
			fi:
		    od:
		od:

		# create the legal_positions matrix
		legal_positions:=linalg[matrix](2,linalg[coldim](ambientmatrix)):
		for aa from nops(lambda)+1 to linalg[coldim](ambientmatrix) do
		    lambda:=[op(lambda),0]:
		od:
		for aa from nops(nu)+1 to linalg[coldim](ambientmatrix) do
		    nu:=[op(nu),0]:
		od:

		lower_boundary:=[]:
		for aa from 1 to linalg[coldim](ambientmatrix) do
		    lower_boundary:=[op(lower_boundary),aa-1]:
		od:
		
		diff:=nu-lambda:

		for aa from 1 to linalg[coldim](ambientmatrix) do
		    if(op(aa,diff)=0) then
			legal_positions[1,aa]:=infinity:
			legal_positions[2,aa]:=infinity:

		    else
			legal_positions[1,aa]:=op(aa,lower_boundary)+op(aa,nu):
			legal_positions[2,aa]:=legal_positions[1,aa]-op(aa,diff)+1:

			legal_positions[1,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[1,aa]+1:
			legal_positions[2,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[2,aa]+1:
		    fi:
		od:
	fi:

	# odd quadric case
	if(flagvarietytype=oddQ) then
		ambientmatrix:=linalg[matrix](1,2*Lie_index-1):
		for aa from 1 to linalg[coldim](ambientmatrix) do
		    ambientmatrix[1,aa]:=_:
		od:
		# create the legal_positions matrix
		legal_positions:=linalg[matrix](2,linalg[coldim](ambientmatrix)):
		for aa from nops(lambda)+1 to linalg[coldim](ambientmatrix) do
		    lambda:=[op(lambda),0]:
		od:
		for aa from nops(nu)+1 to linalg[coldim](ambientmatrix) do
		    nu:=[op(nu),0]:
		od:
		
		diff:=nu-lambda:

		for aa from 1 to linalg[coldim](ambientmatrix) do
		    if(op(aa,diff)=0) then
			legal_positions[1,aa]:=infinity:
			legal_positions[2,aa]:=infinity:

		    else
			legal_positions[1,aa]:=op(aa,nu):
			legal_positions[2,aa]:=legal_positions[1,aa]-op(aa,diff)+1:
			legal_positions[1,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[1,aa]+1:
			legal_positions[2,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[2,aa]+1:
		    fi:
		od:

	fi:

	# even quadric case
	if(flagvarietytype=evenQ) then
		ambientmatrix:=linalg[matrix](2,2*(Lie_index-2)):
		for aa from 1 to 2 do
		  for bb from 1 to linalg[coldim](ambientmatrix) do
		      if(aa=1) then
		           if(bb>=Lie_index-2) then
		               ambientmatrix[aa,bb]:=_:
                           else
			       ambientmatrix[aa,bb]:=Z:
			   fi:
		      else
		           if(bb<=Lie_index-1) then
		               ambientmatrix[aa,bb]:=_:
                           else
			       ambientmatrix[aa,bb]:=Z:
			   fi:			
		      fi:
		  od:
		od:

		# create the legal_positions matrix
		legal_positions:=linalg[matrix](2,linalg[coldim](ambientmatrix)):
		for aa from nops(lambda)+1 to linalg[coldim](ambientmatrix) do
		    lambda:=[op(lambda),0]:
		od:
		for aa from nops(nu)+1 to linalg[coldim](ambientmatrix) do
		    nu:=[op(nu),0]:
		od:

		lower_boundary:=[]:
		for aa from 1 to linalg[coldim](ambientmatrix) do
		    if(aa<=Lie_index-1) then
         		    lower_boundary:=[op(lower_boundary),0]:
                    else
                            lower_boundary:=[op(lower_boundary),1]:
		    fi:
		od:
		
		diff:=nu-lambda:

		for aa from 1 to linalg[coldim](ambientmatrix) do
		    if(op(aa,diff)=0) then
			legal_positions[1,aa]:=infinity:
			legal_positions[2,aa]:=infinity:

		    else
			legal_positions[1,aa]:=op(aa,lower_boundary)+op(aa,nu):
			legal_positions[2,aa]:=legal_positions[1,aa]-op(aa,diff)+1:
			legal_positions[1,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[1,aa]+1:
			legal_positions[2,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[2,aa]+1:
		    fi:
		od:
	fi:



	if(flagvarietytype=E and Lie_index=6) then
		ambientmatrix:=linalg[matrix](4,8):
		for aa from 1 to 4 do
		    for bb from 1 to 8 do
			if(aa=1) then
				if(bb>=4 and bb<=8) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
			if(aa=2) then
				if(bb>=4 and bb<=6) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
			if(aa=3) then
				if(bb>=3 and bb<=5) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
			if(aa=4) then
				if(bb>=1 and bb<=5) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
		    od:
		od:
		
		# create the legal_positions matrix
		legal_positions:=linalg[matrix](2,8):
		lower_boundary:=[0,0,0,0,0,2,3,3]:
		for aa from nops(lambda)+1 to 8 do
		    lambda:=[op(lambda),0]:
		od:
		for aa from nops(nu)+1 to 8 do
		    nu:=[op(nu),0]:
		od:
		
		diff:=nu-lambda:

		for aa from 1 to 8 do
		    if(op(aa,diff)=0) then
			legal_positions[1,aa]:=infinity:
			legal_positions[2,aa]:=infinity:

		    else
			legal_positions[1,aa]:=op(aa,lower_boundary)+op(aa,nu):
			legal_positions[2,aa]:=legal_positions[1,aa]-op(aa,diff)+1:

			legal_positions[1,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[1,aa]+1:
			legal_positions[2,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[2,aa]+1:
		    fi:
		od:
	fi:

	if(flagvarietytype=E and Lie_index=7) then
		ambientmatrix:=linalg[matrix](9,9):
		for aa from 1 to 9 do
		    for bb from 1 to 9 do
			if(aa=1) then
				if(bb>=9 and bb<=9) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
			if(aa=2) then
				if(bb>=9 and bb<=9) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
			if(aa=3) then
				if(bb>=9 and bb<=9) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
			if(aa=4) then
				if(bb>=8 and bb<=9) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
			if(aa=5) then
				if(bb>=5 and bb<=9) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
			if(aa=6) then
				if(bb>=5 and bb<=9) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
			if(aa=7) then
				if(bb>=5 and bb<=7) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
			if(aa=8) then
				if(bb>=4 and bb<=6) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
			if(aa=9) then
				if(bb>=1 and bb<=6) then
				   ambientmatrix[aa,bb]:=_:
				else
				   ambientmatrix[aa,bb]:=Z:
				fi:
			fi:
		    od:
		od:
		
		# create the legal_positions matrix
		legal_positions:=linalg[matrix](2,9):
		lower_boundary:=[0,0,0,0,0,0,2,3,3]:
		for aa from nops(lambda)+1 to 9 do
		    lambda:=[op(lambda),0]:
		od:
		for aa from nops(nu)+1 to 9 do
		    nu:=[op(nu),0]:
		od:
		
		diff:=nu-lambda:

		for aa from 1 to 9 do
		    if(op(aa,diff)=0) then
			legal_positions[1,aa]:=infinity:
			legal_positions[2,aa]:=infinity:

		    else
			legal_positions[1,aa]:=op(aa,lower_boundary)+op(aa,nu):
			legal_positions[2,aa]:=legal_positions[1,aa]-op(aa,diff)+1:

			legal_positions[1,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[1,aa]+1:
			legal_positions[2,aa]:=linalg[rowdim](ambientmatrix)-legal_positions[2,aa]+1:
		    fi:
		od:
	fi:
	
	# recursively create all SYT of shape \nu/\lambda
	# and then check the theorem
	ans:=createSYT(mu,ambientmatrix,legal_positions,maxnumberinserted+1, printoutput):

	# returns the answer and the rectified tableaux

	# now tack on the 2^(shortroots statistic)

	if(flagvarietytype=LG or flagvarietytype=oddQ) then
		ans:=ans*2^(shortroots(nu,flagvarietytype,Lie_index)-shortroots(lambda,flagvarietytype,Lie_index)-shortroots(mu,flagvarietytype,Lie_index)):
	fi:
	return(ans):
end:

# recursively create all SYT of shape \nu/\lambda. Once a SYT is
# fully created, we apply jdt to compute the result

`createSYT`:=proc(mu, ambientmatrixin, legalpositionsin,lastnumberinserted, printoutput)
	local ii, ambientmatrix, ans, legalpositions,
              tempambientmatrix;

	ans:=0:
	ambientmatrix:=linalg[matrix](ambientmatrixin):
	legalpositions:=linalg[matrix](legalpositionsin):

	if(lastnumberinserted>1) then
		# idea is to go through all columns and look for
		# a _corner_ to put the current number, and then branch
		for ii from linalg[coldim](legalpositions) to 1 by -1 do
	            ambientmatrix:=linalg[matrix](ambientmatrixin):
	            legalpositions:=linalg[matrix](legalpositionsin):		    
		    if(legalpositions[1,ii]<infinity) then
		      # check if we have a corner (i.e., we're in column 1,
		      # or the previous column is shorter)
		      if(ii=linalg[coldim](legalpositions)) then
			 ambientmatrix[legalpositions[1,ii],ii]:=lastnumberinserted-1:
			 if(legalpositions[2,ii]=legalpositions[1,ii]) then
				legalpositions[2,ii]:=infinity:
				legalpositions[1,ii]:=infinity:
			 else
				legalpositions[1,ii]:=legalpositions[1,ii]+1:			 
			 fi:
			 ans:=ans+createSYT(mu, ambientmatrix, legalpositions, 
                                                lastnumberinserted-1, printoutput):
		      else if(legalpositions[1,ii+1]>legalpositions[1,ii]) then
			      ambientmatrix[legalpositions[1,ii],ii]:=lastnumberinserted-1:
			      if(legalpositions[2,ii]=legalpositions[1,ii]) then
				legalpositions[2,ii]:=infinity:
				legalpositions[1,ii]:=infinity:
			      else
				legalpositions[1,ii]:=legalpositions[1,ii]+1:			 
			      fi:
			      ans:=ans+createSYT(mu, ambientmatrix, legalpositions, 
                                                     lastnumberinserted-1, printoutput):
                           fi:
		      fi:
		    fi:

		od:
	else # apply jdt and then verify if it counts
	        tempambientmatrix:=linalg[matrix](ambientmatrix):
	        ambientmatrix:=cominjdt(ambientmatrix):
		ans:=LRcheck(mu,ambientmatrix):

		if(ans=1 and printoutput=1) then
			  print(tempambientmatrix, ambientmatrix);
		fi:
	fi:
	return(ans);
end:


`cominjdt`:=proc(ambientmatrixin)
	local ii,jj, ambientmatrix, flag, maxboxfound,slideposition,
              jdtslidedone, rightvalue, upvalue, seen_hyphen_flag, endrow, endcol:

	ambientmatrix:=linalg[matrix](ambientmatrixin):

	maxboxfound:=1:
	while(maxboxfound=1) do
		maxboxfound:=0:
		# the idea is that we look for all pairs of "dominos" consisting of a 
		# "-" below or to the right of a number. Store the current one, and 
		# if there's a bigger one, replace
		
		slideposition:=[infinity,-infinity]:
		for ii from 1 to linalg[rowdim](ambientmatrix) do
		    for jj from 1 to linalg[coldim](ambientmatrix) do
			if(ambientmatrix[ii,jj]=_) then
				if(ii>1) then
					 if(ambientmatrix[ii-1,jj]<>_ and ambientmatrix[ii-1,jj]<>Z) then
						maxboxfound:=1: # may not be a max box, but one exists
						if(op(1,slideposition)>=ii and op(2,slideposition)<=jj) then
							slideposition:=[ii,jj]:
						fi:
					 fi:
				fi:
				if(jj<linalg[coldim](ambientmatrix)) then
					 if(ambientmatrix[ii,jj+1]<>_ and ambientmatrix[ii,jj+1]<>Z) then
						maxboxfound:=1:
						if(op(1,slideposition)>=ii and op(2,slideposition)<=jj) then
							slideposition:=[ii,jj]:
						fi:						
					 fi:
				fi:
			fi:
		    od:
		od:


	        # at this point, we've found a maximal box underneath \nu/\lambda to slide into
		
		if(maxboxfound=1) then
		      jdtslidedone:=0:
		      while(jdtslidedone=0) do
		         if((op(2,slideposition)<linalg[coldim](ambientmatrix)) and 
                            (ambientmatrix[op(1,slideposition),op(2,slideposition)+1]<>Z) and
			    (ambientmatrix[op(1,slideposition),op(2,slideposition)+1]<>_)) then				                rightvalue:=ambientmatrix[op(1,slideposition),op(2,slideposition)+1]: 
		         else
			    rightvalue:=infinity:
		         fi:

		         if((op(1,slideposition)>1) and 
                            (ambientmatrix[op(1,slideposition)-1,op(2,slideposition)]<>Z) and
			    (ambientmatrix[op(1,slideposition)-1,op(2,slideposition)]<>_)) then				                upvalue:=ambientmatrix[op(1,slideposition)-1,op(2,slideposition)]: 
		         else
			    upvalue:=infinity:
		         fi:

			 if((rightvalue=infinity) and (upvalue=infinity)) then
			    jdtslidedone:=1:

			 else # now do the slide
			    if(upvalue<rightvalue) then
				ambientmatrix[op(1,slideposition),op(2,slideposition)]:=
				      upvalue:
				ambientmatrix[op(1,slideposition)-1,op(2,slideposition)]:=_:
				slideposition:=[op(1,slideposition)-1,op(2,slideposition)]:
			    else
				ambientmatrix[op(1,slideposition),op(2,slideposition)]:=
				      rightvalue:
				ambientmatrix[op(1,slideposition),op(2,slideposition)+1]:=_:
				slideposition:=[op(1,slideposition),op(2,slideposition)+1]:
			    fi:
			 fi:
		      od:
		fi:
	od:
	return(ambientmatrix):
end:

# given ambient matrix and mu, we check if the filling is of straight shape mu AND
# the consecutive tableau that has 1,2,3,...i in the first column, i+1,i+2,... in the
# second etc
 
`LRcheck`:=proc(muin, ambientmatrixin)
	local ii, jj, ambientmatrix, failure, mu, 
              currentlabel, numberincolumn, colended:

	ambientmatrix:=linalg[matrix](ambientmatrixin):
	ii:=linalg[rowdim](ambientmatrix):
	jj:=1:
	failure:=0:
	mu:=muin:

	
	for ii from nops(mu)+1 to linalg[coldim](ambientmatrix) do
	    mu:=[op(mu),0]:
	od:

	jj:=1:
	currentlabel:=1:
	while(jj<=linalg[coldim](ambientmatrix) and failure=0) do
		# Check that the number of numbers in each column matches that of mu
		# check that the numbers are in consecutive order going up columns	
		numberincolumn:=0:
		colended:=0:
		ii:=linalg[rowdim](ambientmatrix):
		while(colended=0 and ii>=1 and failure=0) do
			if((ambientmatrix[ii,jj]<>Z) and (ambientmatrix[ii,jj]<>_)) then
				if(ambientmatrix[ii,jj]=currentlabel) then
					currentlabel:=currentlabel+1:
					numberincolumn:=numberincolumn+1:
					if(ii>1) then
						 if((ambientmatrix[ii-1,jj]=Z) or (ambientmatrix[ii-1,jj]=_)) then
							colended:=1:
						 fi:
					fi:
				else
				        failure:=1:
				fi:
			fi: 
			ii:=ii-1:	    
		od: 	
		if(numberincolumn<>op(jj, mu)) then
			failure:=1:
		fi:
		jj:=jj+1:
	od:

	if(failure=1) then
		return(0):
	else
	        return(1):
	fi:
end:

`shortroots`:=proc(theshapein, flagvarietytype, Lie_index)
	local ans, ii,jj, theshape:

	ans:=0:

	theshape:=theshapein:

	if(flagvarietytype=oddQ) then
		for ii from nops(theshape)+1 to 2*Lie_index-1 do
		    theshape:=[op(theshape),0]:
		od:

		if(op(Lie_index,theshape)=1) then
			ans:=1:
		fi:
	fi:	
	
	if(flagvarietytype=LG) then

		for ii from nops(theshape)+1 to Lie_index do
		    theshape:=[op(theshape),0]:
		od:

		for jj from 1 to Lie_index do
		    ans:=ans+max(0,op(jj,theshape)-1):
		od:
			       
	fi:
	return(ans);
end: