print(`This is TASEP: A Maple package for calculating various quantities`):
print(`related to Totally Asymmetric Simple Exclusion Processes.`):
print(`Written by Arvind Ayyer, Rutgers University, October 20 2007.`):
print(`These include open systems on a one-dimensional lattice, closed systems`):
print(`on a ring, and even restricted systems on a one-dimensional lattice.`):
print(`For a list of procedures, type Help();`):
print(`For a list of supplementary procedures, type Help1();`):




Help := proc():

if args = NULL then

	print(`TASEP contains the following main procedures:`):
	print(``):
	print(`For steady state probabilities: SteadyState, ClosedSteadyState, RestrictedSteadyState,`):
	print(`For densities: Opavoc, Clavoc, Ravoc`):
	print(`For two point correlations: Opcorr, Clocorr, Recorr`):

fi:

if nops([args]) = 1 and op(1,[args]) = `SteadyState` then

	print(`SteadyState(sites,colours,Tr1,Tr2,Tr3,p) calculates`):
	print(`the probabilities of each state when the TASEP with`):
	print(`#positions = sites, #particles = colours`):
	print(`is in its steady state where Tr1 is the vector rate of entering,`):
	print(`Tr2 is the matrix rate in the bulk and`):
	print(`Tr3 is the vector rate of exiting`):
	print(`p is a symbolic variable for the probabilities.`):
	print(`For example, try SteadyState(2,2,GenMat(2,[x]),GenMat(2,[1]),GenMat(2,[y]),p);`):
fi:


if nops([args]) = 1 and op(1,[args]) = `ClosedSteadyState` then

	print(`ClosedSteadyState(colours,config,Tr,p) calculates`):
	print(`the steady state probabilities of each state for the closed  TASEP where`):
	print(`config is the list of number of particles of each color`):
	print(`Tr is the matrix rate and`):
	print(`p is a symbolic variable for the probabilities.`):
	print(`For example, try ClosedSteadyState(3,[1,1,1],GenMat2(3,[1$2]),p);`):
fi:



if nops([args]) = 1 and op(1,[args]) = `RestrictedSteadyState` then

	print(`RestrictedSteadyState(sites,colours,midcols,x,Tr2,y,k)`):
	print(`calculates the steady state probabilities of the Restricted TASEP`):
	print(`where x is the rate of entering for colours-1 over 0,`):
	print(`Tr2 is the matrix rate in the bulk and Tr3 is the rate`):
	print(`of exiting for colours-1 leaving 0, and midcols is the`):
	print(`list consisting of the number of particles of type `):
	print(`colours 1,...,colours-2 which are fixed in the system.`):
	print(`For example, try RestrictedSteadyState(2,3,[1],x,GenMat2(3,[1,1]),y,p);`):

fi:

if nops([args]) = 1 and op(1,[args]) = `Opavoc` then

	print(`Opavoc(sites,colours,alpha,beta,gamma,c,i) calculates`):
	print(`the average occupation for colour c on site i for the open TASEP `):
	print(`where alpha, beta and gamma are the rates.`):
	print(`For example, try`):
	print(`seq(Opavoc(2,2,GenMat2(2,[x$1]),GenMat2(2,[1$1]),GenMat2(2,[y$1]),1,i),i=1..2);`):


fi:

if nops([args]) = 1 and op(1,[args]) = `Clavoc` then

	print(`Clavoc(colours,config,Tr,c,i) calculates`):
	print(`the average occupation for colour c on site i for the closed TASEP `):
	print(`where config and Tr are as in ClosedSteadyState `):
	print(`For example, try`):
	print(`seq(Clavoc(3,[1,1,1],GenMat2(3,[1$2]),2,i),i=1..3);`):


fi:


if nops([args]) = 1 and op(1,[args]) = `Ravoc` then

	print(`Ravoc(sites,colours,midcols,x,Tr2,y,c,i) calculates`):
	print(`the average density in the steady state of colour c on site i`):
	print(`for the  Restricted TASEP process`):
	print(`where the quantities are as in RestrictedSteadyState`):
	print(`For example, try`):
	print(`seq(Ravoc(3,3,[1],x,GenMat2(3,[1$2]),y,2,i),i=1..3);`):


fi:

if nops([args]) = 1 and op(1,[args]) = `Opcorr` then

	print(`Opcorr(sites,colours,midcols,x,Tr2,y,c1,n1,c2,n2) calculates`):
	print(`the correlation in the steady state between colour c1 on site n1`):
	print(`and colour c2 on site n2 for the open TASEP process`):
	print(`where the quantities are as in SteadyState.`):
	print(`For example, try`):
	print(`Opcorr(2,3,GenMat2(3,[x$2]),GenMat2(3,[1$2]),GenMat2(3,[y$2]),0,1,0,2);`):

fi:
if nops([args]) = 1 and op(1,[args]) = `Clocorr` then

	print(`Clocorr(colours,config,Tr,c1,n1,c2,n2) calculates`):
	print(`the correlation in the steady state between colour c1 on site n1`):
	print(`and colour c2 on site n2 for the Closed TASEP process`):
	print(`where the quantities are as in ClosedSteadyState.`):
	print(`For example, try`):
	print(`Clocorr(3,[1,1,1],GenMat2(3,[1$2]),0,1,1,3);`):

fi:


if nops([args]) = 1 and op(1,[args]) = `Recorr` then

	print(`Recorr(sites,colours,midcols,x,Tr2,y,c1,n1,c2,n2) calculates`):
	print(`the correlation in the steady state between colour c1 on site n1`):
	print(`and colour c2 on site n2 for the  Restricted TASEP process`):
	print(`where the quantities are as in RestrictedSteadyState.`):
	print(`For example, try`):
	print(`Recorr(3,3,[1],x,GenMat2(3,[1$2]),y,0,1,0,2);`):

fi:


end:

Help1 := proc():

if args = NULL then

	print(`TASEP contains the following supplementary procedures:`):
	print(``):
	print(`Graphs for the processes`):
	print(`OpenGraph, ClosedGraph, RestrictedOpenGraph`):
	print(`Matrices for defining boundary rates`):
	print(`GenMat, GenMat2, GenMat3 `):
	print(`Catalan triangle numbers: cattri`):

fi:

if nops([args]) = 1 and op(1,[args]) = `OpenGraph` then

	print(`OpenGraph(n,c) creates the relevant graph for the open TASEP`):
	print(`with n sites and c colours (labelled 0,..,c-1) assuming the higher colour`):
	print(`can be exchanged with the lower colour to its right in the bulk.`):
	print(`It returns the set of vertices followed by the set of edges.`):
	print(`For example, try OpenGraph(2,2);`):
fi:

if nops([args]) = 1 and op(1,[args]) = `ClosedGraph` then
	
	print(`ClosedGraph(n,c) creates the relevant graph for the closed TASEP`):
	print(`with n sites and c colours (labelled 0,..,c-1) assuming the higher colour`):
	print(`can be exchanged with the lower colour to its right.`):
	print(`It returns the set of vertices followed by the set of edges.`):
	print(`For example, try ClosedGraph(2,2);`):

fi:

if nops([args]) = 1 and op(1,[args]) = `RestrictedOpenGraph` then

	print(`RestrictedOpenGraph(n,c,midcols) creates the relevant graph for the restricted TASEP`):
	print(`with n sites and c colours (labelled 0,..,c-1) assuming the higher colour`):
	print(`can be exchanged with the lower colour to its right in the bulk, but only the highest`):
	print(`and the lowest color can move out of the boundary and the list midcols of size c-2`):
	print(`contains the numbers of the mid colours in the system which are fixed.`):
	print(`It returns the set of vertices followed by the set of edges.`):
	print(`For example, try RestrictedOpenGraph(2,3,[1]);`):

fi:

if nops([args]) = 1 and op(1,[args]) = `GenMat` then

	print(`GenMat(c,A) given an integer c (the number of colours) and a list A of size c-1`):
	print(`returns a matrix of rates where colour i exchanges with colour j by the rate`):
	print(`given by sum(A[k],k=i..j-1). For example, try GenMat(3,[a,b]);`):

fi:

if nops([args]) = 1 and op(1,[args]) = `GenMat2` then

	print(`GenMat2(c,A) given an integer c (the number of colours) and a list A of size c-1`):
	print(`returns a matrix of rates where colour i exchanges with colour j by the rate`):
	print(`given by A[j-i]. For example, try GenMat2(3,[a,b]);`):

fi:

if nops([args]) = 1 and op(1,[args]) = `GenMat3` then

	print(`GenMat3(c,A) given an integer c (the number of colours) and a list of lists A `):
	print(`of size (c-1) times (c-1) returns a matrix of rates where colour i exchanges with`):
	print(`colour j by the rate given by A[i][j]. For example try GenMat3(3,[[0,a,b],[a,0,c],[b,c,0]]);`):
	
fi:

if nops([args]) = 1 and op(1,[args]) = `cattri` then

	print(`cattri(n,m) returns the catalan triangle number C^n_m defined for m=0..n`):
	print(`and n>=0 which is given by the expression binomial(m+n,n)*(n-m+1)/(n+1).`):
	print(`cattri(n,n) returns the nth catalan number. For example, try cattri(4,2);`):

fi:

end:






#################################################GRAPH CREATION############################################################


#OpenGraph(sites,colours) creates the relevant graph for the open system with edges(0,...,n-1)^d and vertices for transition elements assuming the higher colour can be exchanged with the lower colour in the bulk as well as the boundary
OpenGraph := proc(sites,colours) local i,j,k,Ve2,Ed,Ve,tem:

Ve2 := { seq(convert(i,base,colours) ,i=0..colours^sites-1) }:

Ve := {}:

for i from 1 to nops(Ve2) do

	if nops(Ve2[i]) < sites then
		Ve := Ve union {[op(Ve2[i]),0$(sites-nops(Ve2[i]))]}:
	else
		Ve := Ve union {Ve2[i]}:
	fi:
od:

Ed := {}:

for i from 1 to colours^sites do

	tem := Ve[i][1]:
	Ed := Ed union {seq([Ve[i],[k,op(2..sites,Ve[i])]],k=tem+1..colours-1)}:

	tem := Ve[i][sites]:
	Ed := Ed union {seq([Ve[i],[op(1..sites-1,Ve[i]),k]],k=0..tem-1) }:


	for j from 1 to sites-1 do

		if Ve[i][j] > Ve[i][j+1] then
			Ed := Ed union { [Ve[i], [op(1..j-1,Ve[i]), Ve[i][j+1], Ve[i][j], op(j+2..sites,Ve[i])]] }:
		fi:

	od:
od:

return Ve,Ed:

end:
	


#ClosedGraph(sites,colours) creates the relevant graph for the closed system with edges(0,...,n-1)^d and vertices for transition elements assuming the higher colour can be exchanged with the lower colour in the bulk (there is no boundary)
ClosedGraph := proc(sites,colours) local i,j,k,Ve2,Ed,Ve,tem:

Ve2 := { seq(convert(i,base,colours) ,i=0..colours^sites-1) }:

Ve := {}:

for i from 1 to nops(Ve2) do

	if nops(Ve2[i]) < sites then
		Ve := Ve union {[op(Ve2[i]),0$(sites-nops(Ve2[i]))]}:
	else
		Ve := Ve union {Ve2[i]}:
	fi:
od:

Ed := {}:

for i from 1 to colours^sites do

	if Ve[i][1] < Ve[i][sites] then
		Ed := Ed union {[Ve[i],[Ve[i][sites],op(2..sites-1,Ve[i]),Ve[i][1]]] }:
	fi:


	for j from 1 to sites-1 do

		if Ve[i][j] > Ve[i][j+1] then
			Ed := Ed union { [Ve[i], [op(1..j-1,Ve[i]), Ve[i][j+1], Ve[i][j], op(j+2..sites,Ve[i])]] }:
		fi:

	od:
od:

return Ve,Ed:

end:

#RestrictedOpenGraph(sites,colours,midcols) creates the relevant graph for the open system with edges(0,...,n-1)^d and vertices for transition elements assuming the higher colour can be exchanged with the lower colour in the bulk, but only the highest and the lowest color can move out of the boundary and the list midcols of size k-2 contains the numbers of the mid colours in the system
RestrictedOpenGraph := proc(sites,colours,midcols) local i,j,k,Ve2,Ed,Ve,tem,collist:

if sites < 2 then
	ERROR(`Need atleast two sites`):
fi:

if add(midcols[i],i=1..nops(midcols)) > sites then
	ERROR(`Too few sites to accommodate all the middle colours`):
fi:

Ve2 := { seq(convert(i,base,colours) ,i=0..colours^sites-1) }:

Ve := {}:

for i from 1 to nops(Ve2) do

	if nops(Ve2[i]) < sites then
		Ve := Ve union {[op(Ve2[i]),0$(sites-nops(Ve2[i]))]}:
	else
		Ve := Ve union {Ve2[i]}:
	fi:
od:

Ve2 := Ve:
Ve := {}:


for i from 1 to nops(Ve2) do
	collist := [0$(colours-2)]:
	for j from 1 to sites do
		if Ve2[i][j] <>0 and Ve2[i][j] <> colours-1 then
			collist[Ve2[i][j]]:=collist[Ve2[i][j]]+1:
		fi:
	od:

	if collist = midcols then
		Ve := Ve union {Ve2[i]}:
	fi:
od:
	
Ed := {}:

for i from 1 to nops(Ve) do

	if Ve[i][1] = 0 then
		Ed := Ed union {[Ve[i],[colours-1,op(2..sites,Ve[i])]]}:
	fi:

	if Ve[i][sites] = colours-1 then
		Ed := Ed union {[Ve[i],[op(1..sites-1,Ve[i]),0]]}:
	fi:


	for j from 1 to sites-1 do

		if Ve[i][j] > Ve[i][j+1] then
			Ed := Ed union { [Ve[i], [op(1..j-1,Ve[i]), Ve[i][j+1], Ve[i][j], op(j+2..sites,Ve[i])]] }:
		fi:

	od:
od:

return Ve,Ed:

end:
	


#################################################STEADY STATE PROBABILITIES############################################################

#SimpleSteadyState(sites,colours,alpha,beta,gamma,k) calculates the probabilities of each state when the TASEP is in its steady state where alpha is the rate of entering, beta is the rate in the bulk and gamma is the rate of exiting.
SimpleSteadyState := proc(sites, colours, alpha, beta, gamma,k) local i,j,k2,Gr,Ed,Ve,eqs,temp:

Gr := OpenGraph(sites,colours):

Ed := Gr[2]:
Ve := Gr[1]:

k2 := {}:

for i from 1 to nops(Ve) do
	k2 := k2 union {k[op(Ve[i])]}:
od:

eqs := {}:

for i from 1 to nops(Ve) do
	temp := 0:
	for j from 1 to nops(Ed) do
		if Ve[i] = Ed[j][1] then
			if sites = 1 then
				if op(Ed[j][1]) < op(Ed[j][2]) then
					temp := temp - alpha*k[op(Ed[j][1])]:
				else
					temp := temp - gamma*k[op(Ed[j][1])]:
				fi:
			else
				if Ed[j][1][1] <> Ed[j][2][1] and Ed[j][1][2] = Ed[j][2][2] then
					temp := temp - alpha*k[op(Ed[j][1])]:
				elif Ed[j][1][sites] <> Ed[j][2][sites] and Ed[j][1][sites-1] = Ed[j][2][sites-1] then
					temp := temp - gamma*k[op(Ed[j][1])]:
				else
					temp := temp - beta*k[op(Ed[j][1])]:
				fi:
			fi:
		fi:

		if Ve[i] = Ed[j][2] then
			
			if sites = 1 then
				if op(Ed[j][1]) < op(Ed[j][2]) then
					temp := temp + alpha*k[op(Ed[j][1])]:
				else
					temp := temp + gamma*k[op(Ed[j][1])]:
				fi:
			else
				if Ed[j][1][1] <> Ed[j][2][1] and Ed[j][1][2] = Ed[j][2][2] then
					temp := temp + alpha*k[op(Ed[j][1])]:
				elif Ed[j][1][sites] <> Ed[j][2][sites] and Ed[j][1][sites-1] = Ed[j][2][sites-1] then
					temp := temp + gamma*k[op(Ed[j][1])]:
				else
					temp := temp + beta*k[op(Ed[j][1])]:
				fi:
			fi:
		fi:
	od:
	eqs := eqs union {temp}:
od:

eqs := eqs union {add(k2[i],i=1..nops(k2))-1}:


return solve(eqs,k2):



end:


#SteadyState(sites,colours,Tr1,Tr2,Tr3,k) calculates the probabilities of each state when the TASEP is in its steady state where Tr1 is the matrix rate of entering, Tr2 is the matrix rate in the bulk and Tr3 is the matrix rate of exiting.
#Tr1, Tr2 and Tr3 must be list of lists
SteadyState := proc(sites, colours, Tr1, Tr2, Tr3, k) local i,j,l,k2,Gr,Ed,Ve,eqs,temp:

Gr := OpenGraph(sites,colours):

Ed := Gr[2]:
Ve := Gr[1]:

k2 := {}:

for i from 1 to nops(Ve) do
	k2 := k2 union {k[op(Ve[i])]}:
od:

eqs := {}:

for i from 1 to nops(Ve) do
	temp := 0:
	for j from 1 to nops(Ed) do
		if Ve[i] = Ed[j][1] then
			if sites = 1 then
				if op(Ed[j][1]) < op(Ed[j][2]) then
					temp := temp - Tr1[Ed[j][2][1]+1][Ed[j][1][1]+1]*k[op(Ed[j][1])]:
				else
					temp := temp - Tr3[Ed[j][1][1]+1][Ed[j][2][1]+1]*k[op(Ed[j][1])]:
				fi:
			else
				if Ed[j][1][1] <> Ed[j][2][1] and Ed[j][1][2] = Ed[j][2][2] then
					temp := temp - Tr1[Ed[j][2][1]+1][Ed[j][1][1]+1]*k[op(Ed[j][1])]:
				elif Ed[j][1][sites] <> Ed[j][2][sites] and Ed[j][1][sites-1] = Ed[j][2][sites-1] then
					temp := temp - Tr3[Ed[j][1][sites]+1][Ed[j][2][sites]+1]*k[op(Ed[j][1])]:
				else
					for l from 1 to sites-1 do
						if Ed[j][1][l] <> Ed[j][2][l] and Ed[j][1][l+1] <> Ed[j][2][l+1] then
							temp := temp - Tr2[Ed[j][2][l]+1][Ed[j][1][l]+1]*k[op(Ed[j][1])]:
						fi:
					od:	
				fi:
			fi:
		fi:

		if Ve[i] = Ed[j][2] then
			if sites = 1 then
				if op(Ed[j][1]) < op(Ed[j][2]) then
					temp := temp + Tr1[Ed[j][2][1]+1][Ed[j][1][1]+1]*k[op(Ed[j][1])]:
				else
					temp := temp + Tr3[Ed[j][1][1]+1][Ed[j][2][1]+1]*k[op(Ed[j][1])]:
				fi:
			else

				if Ed[j][1][1] <> Ed[j][2][1] and Ed[j][1][2] = Ed[j][2][2] then
					temp := temp + Tr1[Ed[j][2][1]+1][Ed[j][1][1]+1]*k[op(Ed[j][1])]:
				elif Ed[j][1][sites] <> Ed[j][2][sites] and Ed[j][1][sites-1] = Ed[j][2][sites-1] then
					temp := temp + Tr3[Ed[j][1][sites]+1][Ed[j][2][sites]+1]*k[op(Ed[j][1])]:
				else
					for l from 1 to sites-1 do
						if Ed[j][1][l] <> Ed[j][2][l] and Ed[j][1][l+1] <> Ed[j][2][l+1] then
							temp := temp + Tr2[Ed[j][2][l]+1][Ed[j][1][l]+1]*k[op(Ed[j][1])]:
						fi:
					od:	
				fi:
			fi:
		fi:
	od:
	eqs := eqs union {temp}:
od:

eqs := eqs union {add(k2[i],i=1..nops(k2))-1}:


return solve(eqs,k2):
#return eqs,k2:


end:



#ClosedSteadyState(colours,config,Tr2,k) calculates the probabilities of each state when the TASEP is in its steady state where Tr2 is the matrix rate  and config is a list of length colours with its elements specifying the number of each of colour.
#Tr1, Tr2 and Tr3 must be list of lists
ClosedSteadyState := proc(colours, config, Tr2, k) local i,j,l,k2,Gr,Ed,Ve,Ve2,eqs,temp,sites,cntr:

sites := add(config[i],i=1..nops(config)):

Gr := ClosedGraph(sites,colours):

Ed := Gr[2]:
Ve := {}:
Ve2 := Gr[1]:

k2 := {}:

for i from 1 to nops(Ve2) do
	cntr := [0$colours]:
	for j from 1 to sites do
		cntr[Ve2[i][j]+1] := cntr[Ve2[i][j]+1]+1:
	od:

	if cntr = config then
		Ve := Ve union {Ve2[i]}:
	fi:
od:

for i from 1 to nops(Ve) do
	k2 := k2 union {k[op(Ve[i])]}:
od:

eqs := {}:

for i from 1 to nops(Ve) do
	temp := 0:
	for j from 1 to nops(Ed) do
		if Ve[i] = Ed[j][1] then
			if Ed[j][1][sites] <> Ed[j][2][sites] and Ed[j][1][1] <> Ed[j][2][1] then
				temp := temp - Tr2[Ed[j][2][sites]+1][Ed[j][1][sites]+1]*k[op(Ed[j][1])]:
			fi:

			for l from 1 to sites-1 do
				if Ed[j][1][l] <> Ed[j][2][l] and Ed[j][1][l+1] <> Ed[j][2][l+1] then
					temp := temp - Tr2[Ed[j][2][l]+1][Ed[j][1][l]+1]*k[op(Ed[j][1])]:
				fi:
			od:	

		fi:

		if Ve[i] = Ed[j][2] then
			if Ed[j][1][sites] <> Ed[j][2][sites] and Ed[j][1][1] <> Ed[j][2][1] then
				temp := temp + Tr2[Ed[j][2][sites]+1][Ed[j][1][sites]+1]*k[op(Ed[j][1])]:
			fi:

			for l from 1 to sites-1 do
				if Ed[j][1][l] <> Ed[j][2][l] and Ed[j][1][l+1] <> Ed[j][2][l+1] then
					temp := temp + Tr2[Ed[j][2][l]+1][Ed[j][1][l]+1]*k[op(Ed[j][1])]:
				fi:
			od:	
		fi:
	od:
	eqs := eqs union {temp}:
od:

eqs := eqs union {add(k2[i],i=1..nops(k2))-1}:


return solve(eqs,k2):
#return eqs:


end:


#RestrictedSteadyState(sites,colours,midcols,Tr1,Tr2,Tr3,k) calculates the probabilities of each state when the TASEP is in its steady state where Tr1 is the rate of entering for colours-1 over 0, Tr2 is the matrix rate in the bulk and Tr3 is the rate of exiting for colours-1 leaving 0, and midcols is the list consisting of the number of particles of colours 1,...,colours-2.
#Tr1, Tr2 and Tr3 must be list of lists
#rs := RestrictedSteadyState(6,3,[3],x,GenMat2(4,[1$3]),y,p): {seq(denom(op(2,rs[i])),i=1..nops(rs))};
RestrictedSteadyState := proc(sites, colours, midcols,Tr1, Tr2, Tr3, k) local i,j,l,k2,Gr,Ed,Ve,eqs,temp:

Gr := RestrictedOpenGraph(sites,colours,midcols):

Ed := Gr[2]:
Ve := Gr[1]:

k2 := {}:

for i from 1 to nops(Ve) do
	k2 := k2 union {k[op(Ve[i])]}:
od:

eqs := {}:

for i from 1 to nops(Ve) do
	temp := 0:
	for j from 1 to nops(Ed) do
		if Ve[i] = Ed[j][1] then
			if Ed[j][1][1] <> Ed[j][2][1] and Ed[j][1][2] = Ed[j][2][2] then
				temp := temp - Tr1*k[op(Ed[j][1])]:
			elif Ed[j][1][sites] <> Ed[j][2][sites] and Ed[j][1][sites-1] = Ed[j][2][sites-1] then
				temp := temp - Tr3*k[op(Ed[j][1])]:
			else
				for l from 1 to sites-1 do
					if Ed[j][1][l] <> Ed[j][2][l] and Ed[j][1][l+1] <> Ed[j][2][l+1] then
						temp := temp - Tr2[Ed[j][2][l]+1][Ed[j][1][l]+1]*k[op(Ed[j][1])]:
					fi:
				od:	
			fi:
		fi:

		if Ve[i] = Ed[j][2] then
			if Ed[j][1][1] <> Ed[j][2][1] and Ed[j][1][2] = Ed[j][2][2] then
				temp := temp + Tr1*k[op(Ed[j][1])]:
			elif Ed[j][1][sites] <> Ed[j][2][sites] and Ed[j][1][sites-1] = Ed[j][2][sites-1] then
				temp := temp + Tr3*k[op(Ed[j][1])]:
			else
				for l from 1 to sites-1 do
					if Ed[j][1][l] <> Ed[j][2][l] and Ed[j][1][l+1] <> Ed[j][2][l+1] then
						temp := temp + Tr2[Ed[j][2][l]+1][Ed[j][1][l]+1]*k[op(Ed[j][1])]:
					fi:
				od:	
			fi:
		fi:
	od:
	eqs := eqs union {temp}:
od:

eqs := eqs union {add(k2[i],i=1..nops(k2))-1}:


return solve(eqs,k2):
#return eqs:


end:













#################################################CALCULATIONS BASED ON PROBABILITIES###############################################



#Calculates the common denominator of the stat. prob of the TASEP
Denom := proc(sites,colours,A,B,C) local St,i,p,Den:

St := SteadyState(sites,colours,A,B,C,p):

Den := {seq(denom(op(2,St[i])),i=1..nops(St))}:

return  lcm(op(Den)):

end:






#SimpleAvOccup(sites,colours,alpha,beta,gamma,i) tries to calculate the average occupation on site i using SimpleSteadyState()
SimpleAvOccup := proc(sites,colours,alpha,beta,gamma,i) local j,k,p,St,Av:

St := SimpleSteadyState(sites,colours,alpha,beta,gamma,p):

Av := add(op(2,St[j])*op(i,op(1,St[j])) ,j=1..nops(St)):


return normal(Av):

end:


#AvOccup2(sites,colours,alpha,beta,gamma,i) tries to calculate the average occupation on site i using SteadyState()
#Example: n:=2:c:=2: seq(AvOccup2(n,c,GenMat2(c,[x$(c-1)]),GenMat2(c,[1$(c-1)]),GenMat2(c,[y$(c-1)]),i),i=1..n);
AvOccup2 := proc(sites,colours,A,B,C,i) local j,k,p,St,Av:

St := SteadyState(sites,colours,A,B,C,p):

Av := add(op(2,St[j])*op(i,op(1,St[j])) ,j=1..nops(St)):


return normal(Av):

end:

#Opavoc(sites,colours,alpha,beta,gamma,c,i) tries to calculate the average occupation on site i of the colour c using SteadyState()
#Example: n:=2:c:=2: seq(AvOccup2(n,c,GenMat2(c,[x$(c-1)]),GenMat2(c,[1$(c-1)]),GenMat2(c,[y$(c-1)]),i),i=1..n);
Opavoc := proc(sites,colours,A,B,C,c,i) local j,k,p,St,Av:

if c > colours-1 then
	ERROR(`c has to be less than the number of colours-1`):
fi:

St := SteadyState(sites,colours,A,B,C,p):

Av := 0:

for j from 1 to nops(St) do
	if op(i,op(1,St[j])) = c then
		Av := Av + op(2,St[j]): #Don't multiply by op(i,op(1,St[j]))
	fi:
od:


return normal(Av):

end:


#SimpleNumPart(sites,colours,alpha,beta,gamma,i) calculates the total probability that there are only i particles using SimpleSteadyState() and if particle size = n, it is equivalent to n particles with size = 1
SimpleNumPart := proc(sites,colours,alpha,beta,gamma,i) local j,k,p,St,Num:

St := SimpleSteadyState(sites,colours,alpha,beta,gamma,p):

Num := 0:

for j from 1 to nops(St) do
	if add(op(k,op(1,St[j])),k=1..sites) = i then
		Num := Num + op(2,St[j]):
	fi:
od:

return Num:

end:


#NumPart(sites,colours,A,B,C,i) calculates the total probability that there are only i particles using SteadyState() and if particle size = n, it is equivalent to n particles with size = 1
NumPart := proc(sites,colours,A,B,C,i) local j,k,p,St,Num:

St := SteadyState(sites,colours,A,B,C,p):

Num := 0:

for j from 1 to nops(St) do
	if add(op(k,op(1,St[j])),k=1..sites) = i then
		Num := Num + op(2,St[j]):
	fi:
od:

return Num:

end:

#AvgEta(sites,colours,A,B,C,i,k) calculates the expectation value of the function \eta(i,k) = 1 if colour k is at site i else 0
AvgEta := proc(sites,colours,A,B,C,i,k) local St,p,aveta,j:

St := SteadyState(sites,colours,A,B,C,p):

aveta := 0:

for j from 1 to nops(St) do
	if op(i,op(1,St[j])) = k then
		aveta := aveta + op(2,St[j]):
	fi:
od:

return aveta:

end:


#Corrfn(sites,colours,A,B,C,i,j) calculates the correlation function between i and j
Corrfn := proc(sites,colours,A,B,C,i,j) local p,St,Cf,k:

St := SteadyState(sites,colours,A,B,C,p):

Cf := add(op(2,St[k])*op(i,op(1,St[k]))*op(j,op(1,St[k])) ,k=1..nops(St)):


return normal(Cf):

end:

#Opcorr(sites,colours,A,B,C,c1,n1,c2,n2) calculates the correlation between c1 at n1 and c2 at n2
Opcorr := proc(sites,colours,A,B,C,c1,n1,c2,n2) local p,St,Cf,k:

St := SteadyState(sites,colours,A,B,C,p):

Cf := 0:

for k from 1 to nops(St) do
	if op(n1,op(1,St[k])) = c1 and op(n2,op(1,St[k]))= c2 then
		Cf := Cf + op(2,St[k]):
	fi:
od:

return Cf:

end:





#Current2(sites,colours,A,B,C,i,k) calculates the current of color k at site i using SteadyState()
Current2 := proc(sites,colours,A,B,C,i,k) local cur,St,p,j,m:

if i>sites or k>colours-1 then
	ERROR(`Either site or colour selected is too high`):
fi:

St := SteadyState(sites,colours,A,B,C,p):

cur := 0:


for j from 1 to nops(St) do
	if op(i,op(1,St[j])) = k and op(i+1,op(1,St[j])) < k then
		cur := cur + op(2,St[j]):
	fi:

	if op(i+1,op(1,St[j])) = k and op(i,op(1,St[j])) > k then
		cur := cur - op(2,St[j]):
	fi:

od:

return cur:

end:

#Clavoc(colours,config,Tr2,c,i) tries to calculate the average occupation on site i of the colour c using ClosedSteadyState()
Clavoc := proc(colours,config,Tr2,c,i) local j,k,p,St,Av:

if c > colours-1 then
	ERROR(`c has to be less than the number of colours-1`):
fi:

St := ClosedSteadyState(colours,config,Tr2,p):

Av := 0:

for j from 1 to nops(St) do
	if op(i,op(1,St[j])) = c then
		Av := Av + op(2,St[j]): #Don't multiply by op(i,op(1,St[j]))
	fi:
od:


return normal(Av):

end:


#Clocorr(colours,config,Tr2,c1,n1,c2,n2) calculates the correlation between c1 at n1 and c2 at n2
Clocorr := proc(colours,config,Tr2,c1,n1,c2,n2) local p,St,Cf,k:

St := ClosedSteadyState(colours,config,Tr2,p):

Cf := 0:

for k from 1 to nops(St) do
	if op(n1,op(1,St[k])) = c1 and op(n2,op(1,St[k]))= c2 then
		Cf := Cf + op(2,St[k]):
	fi:
od:

return Cf:

end:





#Ravoc(sites,colours,midcols,x,Tr2,y,c,i) tries to calculate the average occupation on site i of the colour c using RestrictedSteadyState()
Ravoc := proc(sites,colours,midcols,x,Tr2,y,c,i) local j,k,p,St,Av:

if c > colours-1 then
	ERROR(`c has to be less than the number of colours-1`):
fi:

St := RestrictedSteadyState(sites,colours,midcols,x,Tr2,y,p):

Av := 0:

for j from 1 to nops(St) do
	if op(i,op(1,St[j])) = c then
		Av := Av + op(2,St[j]): #Don't multiply by op(i,op(1,St[j]))
	fi:
od:


return normal(Av):

end:


#Recorr(sites,colours,midcols,x,Tr2,y,c1,n1,c2,n2) calculates the correlation function between colour c1 at n1 and colour c2 at n2 for RestrictedSteadyState() divided by the densite of c2 particles at n2
Recorr := proc(sites,colours,midcols,x,Tr2,y,c1,n1,c2,n2) local p,St,Cf,k,i,j:

if c1>colours-1 or c2>colours-1 then
	ERROR(`c1 and c2 have to be within the colours specified`):
fi:

if n1>sites or n2>sites then
	ERROR(`n1 and n2 have to be within the number of sites specified`):
fi:

St := RestrictedSteadyState(sites,colours,midcols,x,Tr2,y,p):

Cf := 0:

for j from 1 to nops(St) do
	if op(n1,op(1,St[j])) = c1 and op(n2,op(1,St[j])) = c2 then
		Cf := Cf + op(2,St[j]):
	fi:
od:

#Cf := Cf/Ravoc(sites,colours,midcols,x,Tr2,y,c2,n2):

return normal(Cf):

end:



#################################################CONJECTURES FOR 2-COLOR SYSTEMS#######################################################





#Denom2col(sites,x,y) calculates the conjectured denominator of the TASEP with sites number of particles
#PROBLEM HERE
Denom2col := proc(sites,x,y) local i,j,den,tem:

den := 0:

for i from 1 to sites do
	tem := 0:
	for j from 1 to sites+2-i do
		tem := tem + x^(sites-j+1)*y^(i+j-2)
	od:

	den := den + cattri(sites-1,i-1)*tem:
od:

return den:

end:







#################################################PLAYING WITH 3-COLOR SYSTEM############################################################
#evaluate the probability of a given configuration given in terms of 0's,1's and 2's in the restricted 3-tasep using the matrix ansatz:
probtau:=proc(L,x,y) local i,j,nums,n,cnt:
option remember:

n:= nops(L):

if n=0 then
	return 1:
fi:

nums := [0$3]:

for i from 1 to n do
	nums[L[i]+1]:=nums[L[i]+1]+1:
od:

cnt:=2*n:

for i from 1 to n-1 do
	if L[i] > L[i+1] and i<cnt then
		cnt := i:
	fi:
od:

#Nondecreasing sequence
if cnt = 2*n then
	return 1/x^nums[1]/y^nums[3]:
fi:

#If not, three cases
if L[cnt]=2 and L[cnt+1]=0 then
	return probtau([op(1..cnt-1,L),2,op(cnt+2..n,L)],x,y)+probtau([op(1..cnt-1,L),0,op(cnt+2..n,L)],x,y):
elif L[cnt]=2 and L[cnt+1]=1 then
	return probtau([op(1..cnt-1,L),1,op(cnt+2..n,L)],x,y):
elif L[cnt]=1 and L[cnt+1]=0 then
	return probtau([op(1..cnt-1,L),1,op(cnt+2..n,L)],x,y):
else
	ERROR(`Something is wrong`):
fi:


end:


#Denom3col(sites,k,x,y) calculates the conjectured denominator of the Restricted 3-TASEP with _sites_ number of particles and _k_ number of 1's.
Denom3col := proc(sites,k,x,y) local i,j,den,tem:

if k > sites or k < 0 then
	return 0:
fi:

den := 0:

if k = 0 then
	for i from 1 to sites do
		tem := 0:
		for j from 1 to sites+2-i do
			tem := tem + x^(sites-j+1)*y^(i+j-2)
		od:

		den := den + cattri(sites-1,i-1)*tem:
	od:
else
	for i from 1 to sites+1-k do
		tem := 0:
		for j from 1 to (sites-k)+2-i do
			tem := tem + x^(sites-k-j+1)*y^(i+j-2) #Total power is sites-k+i-1
		od:

		den := den + cattri(sites+k-1,i-1)*tem:
	od:
fi:

return den:

end:


#recorr3col(sites,k,x,y) calculates the conjectured denominator of the Restricted 3-TASEP with _sites_ number of particles and _k_ number of 1's.
recorr3col := proc(sites,k,x,y) local i,j,den,tem:

if k >= sites or k <= 0 then
	ERROR(`Must have more sites than number of 1's`):
fi:

den := 0:

for i from 1 to sites-k do
	tem := 0:
	for j from 1 to (sites-k)+1-i do
		tem := tem + x^(sites-k-j)*y^(i+j-2) #Total power is sites-k+i-1
	od:
		den := den + cattri(sites+k-1,i-1)*tem:
od:


return x*y*den:

end:


#conjx2(n,n1,i,x,y) returns the conjectured numerator of the 3-tasep of <tau_2>_i
conjx2 := proc(n,n1,i,x,y) local k:

if i <= n1 then
	return factor(add((x*y)^(n-k)*cattri(n-k-1,n-k-1)*Denom3col(k,n1,x,y),k=n1..n-1)):
fi:

if i < n then
	return factor(add((x*y)^k*cattri(k-1,k-1)*Denom3col(n-k,n1,x,y),k=1..n-i) + x^(n-i+1)*Denom3col(i-1,n1,x,y)*add(cattri(n-i-1,k)*y^k,k=0..n-i-1)):
fi:

if i = n then
	return x*Denom3col(n-1,n1,x,y):
fi:

print(`i has to be less than n`):

end:

#sumx2(N,n1,i,x,y) returns the conjectured numerator of the 3-tasep of <tau_2>_i
sumx2 := proc(N,n1,i,x,y) local k,l,s,m1,a,b,mu,m,n,p:

if i <= 1 then
	ERROR(`i must be atleast 2`):
fi:

m:=i-1:
n:=N-m-1:
mu:=n1:

print(`Bulk term`):
print(add(add(add(add(cattri(n+mu-m1+l-s+1,n-mu+m1-b)*cattri(m+m1-1,m-m1-l)*x^(-s)*y^(-b),l=s..m-m1),b=0..n+m1-mu),s=0..m-m1),m1=1..mu-1)):
if m<mu and n<mu then
	print(`Special case`):
#	print(add(add(add(add(cattri(n+mu-m1+l-s+1,n-mu+m1-b)*cattri(m+m1-1,m-m1-l)*x^(-s)*y^(-b),m1=b+mu-n..m-l),l=s..n+m-mu-b),b=0..n+m-mu-s),s=0..n+m-mu)):
	print((x*y)^(N-n1)*add(add(cattri(n+m+mu,n+m-mu-b)*x^(-t)*y^(t-b),t=0..b),b=0..m+n-mu)):
	print(factor(add((x*y)^(N-k)*cattri(N-k-1,N-k-1)*Denom3col(k,mu,x,y),k=mu..N-1))):
fi:


print(`Left boundary term`):
print(add(add(add(cattri(m+mu-1,m-mu-l)*cattri(n+s+1,n-k-1)*x^(s-l)*y^(-k),k=0..n-1),s=0..l),l=0..m-mu) + 
add(add(add(add(cattri(m+mu-1,m-mu-l)*cattri(n+s-t+1,n-k-1)*x^(s-l)*y^(-t-k),k=0..n-1),t=1..s+1),s=0..l),l=0..m-mu) +
add(add(add(cattri(m+mu-1,m-mu-l)*cattri(n-1,n-p)*x^(s-l)*y^(-s-p-1),p=1..n),s=0..l),l=0..m-mu)):



print(`Right boundary term`):
print(add(add(add(add(cattri(m-1,m-p)*cattri(n+mu+p-k,n-mu-a)*x^(-b-k)*y^(b-a),p=k..m),k=0..m),b=0..a),a=0..n-mu)):


return (x*y)^(N-n1)*(

#Bulk term
add(add(add(add(cattri(n+mu-m1+l-s+1,n-mu+m1-b)*cattri(m+m1-1,m-m1-l)*x^(-s)*y^(-b),l=s..m-m1),b=0..n+m1-mu),s=0..m-m1),m1=1..mu-1) + 

#Boundary terms where all 1's are to the left
add(add(add(cattri(m+mu-1,m-mu-l)*cattri(n+s+1,n-k-1)*x^(s-l)*y^(-k),k=0..n-1),s=0..l),l=0..m-mu) + 
add(add(add(add(cattri(m+mu-1,m-mu-l)*cattri(n+s-t+1,n-k-1)*x^(s-l)*y^(-t-k),k=0..n-1),t=1..s+1),s=0..l),l=0..m-mu) + 
add(add(add(cattri(m+mu-1,m-mu-l)*cattri(n-1,n-p)*x^(s-l)*y^(-s-p-1),p=1..n),s=0..l),l=0..m-mu) + 

#Boundary terms where all 1's are to the right
add(add(add(add(cattri(m-1,m-p)*cattri(n+mu+p-k,n-mu-a)*x^(-b-k)*y^(b-a),p=k..m),k=0..m),b=0..a),a=0..n-mu)

):

end:

#This computes the sum from our paper, then multiplies it by the inverse of
#   one of the three different asymptotic forms.  All three scalings are here;
#   two of the three are just commented out with the  #  sign.
asymp := proc(n,rho,x)
   local p,s,y,n1;
   s := 0; n1:=rho*n;
   for p from 0 to n-n1 do 
      y := x^(-p-1)*(p+2*n1)*binomial(2*n-1-p, n+n1-1)/(n+n1);
      s := s+y;
   end do;
   evalf((4/((1+rho)^(1+rho)*(1-rho)^(1-rho)))^(-n)*(2*rho/(1+rho)^2)*(2/(1-rho))*s);      #x=(1-rho)/2    
   #evalf((1/2)*(4/((1+rho)^(1+rho)*(1-rho)^(1-rho)))^(-n)*sqrt(n*Pi*(1-rho^2))*((1)/(rho))*(2*x-1+rho)*s)  #x>(1-rho)/2
   #s*(x*(1-x))^(n+1)*((1-x)/x)^(rho*n)/(1-2*x)           #x<(1-rho)/2
end proc:





#tests a conjecture of gene speer
spconj := proc(m,m1,site) local S,Ve,Ve2,i,j,collist:

Ve2 := { seq(convert(i,base,3) ,i=0..3^m-1) }:
Ve := {}:

for i from 1 to nops(Ve2) do

	if nops(Ve2[i]) < m then
		Ve := Ve union {[op(Ve2[i]),0$(m-nops(Ve2[i]))]}:
	else
		Ve := Ve union {Ve2[i]}:
	fi:
od:

Ve2:=Ve:
Ve:={}:

for i from 1 to nops(Ve2) do
	collist := 0:
	for j from 1 to m do
		if Ve2[i][j] =1 then
			collist := collist+1:
		fi:
	od:

	if collist = m1 and Ve2[i][site]=2 then
		Ve := Ve union {Ve2[i]}:
	fi:
od:

#return Ve:

return add(probtau(Ve[i],x,1),i=1..nops(Ve)):

end:





#tests a conjecture of gene speer and arvind ayyer which is a stronger version of speer's original conjecture
asconj := proc(m,m1,num0,site) local S,Ve,Ve2,i,j,collist,numzero:

Ve2 := { seq(convert(i,base,3) ,i=0..3^m-1) }:
Ve := {}:

for i from 1 to nops(Ve2) do

	if nops(Ve2[i]) < m then
		Ve := Ve union {[op(Ve2[i]),0$(m-nops(Ve2[i]))]}:
	else
		Ve := Ve union {Ve2[i]}:
	fi:
od:

Ve2:=Ve:
Ve:={}:

for i from 1 to nops(Ve2) do
	collist := 0:
	numzero := 0:
	for j from 1 to m do
		if Ve2[i][j] =1 then
			collist := collist+1:
		fi:
		if Ve2[i][j] =0 then
			numzero := numzero+1:
		fi:
	od:

	if collist = m1 and Ve2[i][site]=2 and numzero=num0 then
		Ve := Ve union {Ve2[i]}:
	fi:
od:

#return Ve:

return add(probtau(Ve[i],x,1),i=1..nops(Ve)):

end:



#clrecur(conf,x,y,z) guesses the numerator satisfied by the steady state probability of that configuration of 3 colors
clrecur := proc(conf,x,y,z) local i,gu,nums,conf0,conf2,numzero,numtwo:
option remember:

gu := 0:
nums := [0,0,0]:

for i from 1 to nops(conf) do
	nums[conf[i]+1]:= nums[conf[i]+1]+1:
	
	if i < nops(conf) then
		if conf[i]=2 and conf[i+1]=0 then
			gu := i:
		fi:
	else
		if conf[nops(conf)]=2 and conf[1] =0 then
			gu := nops(conf):
		fi:
	fi:
od:

if nums[2]=0 then
	return 1:
fi:

if nums[1]=0 then
	if nums[3] >= nums[2] then
		return y^(nums[3]-1):
	else
		return 1:
	fi:
fi:

if nums[3]=0 then
	if nums[1] >= nums[2] then
		return y^(nums[1]-1):
	else
		return 1:
	fi:
fi:

if gu = 0 then
	return y^(nums[1]+nums[3]-1):
fi:

if gu < nops(conf) then
	conf0:=[op(1..gu,conf),op(gu+2..nops(conf),conf)]:
	conf2:=[op(1..gu-1,conf),op(gu+1..nops(conf),conf)]:
else
	conf0 := [op(2..nops(conf),conf)]:
	conf2 := [op(1..nops(conf)-1,conf)]:
fi:

return x*clrecur(conf0,x,y,z) + z*clrecur(conf2,x,y,z):


end:






#################################################MATRICES NEEDED TO SPECIFY RATES#######################################################

#GenMat(n,A) generates the "useful" list of lists that could be used in  SteadyState() and other functions depending on a parameter A.
GenMat := proc(n,A) local i,j,t1:

t1 := [[0$n]$n]:


for i from 1 to n do
	for j from i+1 to n do

		t1[j-i][j] := add(A[k],k=j-i..j-1):

	od:
od:

for i from 1 to n do
	for j from 1 to i-1 do
		t1[i][j] := t1[j][i]:
	od:
od:

return t1:

end:


#GenMat2(n,A) generates another "useful" list of lists that could be used in  SteadyState() and other functions depending on a parameter A.
GenMat2 := proc(n,A) local i,j,t1:

t1 := [[0$n]$n]:

for i from 1 to n do
	for j from 1 to n do

		if i <> j then
			t1[i][j] := A[abs(i-j)]:
		fi:

	od:
od:

return t1:

end:


#GenMat3(n,A) generates the simplest list of lists that could be used in  SteadyState() and other functions depending on a list A which is 
GenMat3 := proc(n,A) local i,j,t1:

t1 := [[0$n]$n]:

for i from 1 to n do
	for j from 1 to n do

		if i <> j then
			t1[i][j] := A[i][j]:
		fi:

	od:
od:

return t1:

end:









######################################Generic Useful Programs###################################################
#cattri(n,k) calculates the (n,k)th element of Catalan's triangle where 0<=k<=n
cattri := proc(n,k) option remember:

return binomial(n+k,k)*(n-k+1)/(n+1):

end:


#closfrac(q,denom1) given a fraction q, returns the closest fraction to q whose denominator is not larger than denom1
closfrac := proc(q,denom1) local i,j,d,min1:

d := denom(q):
min1 := numer(q)/d+1:

for i from d+1 to denom1 do
	for j from 1 to i do
		if abs(j/i-q) < abs(min1-q) and j/i <> q then
			min1 := j/i:
		fi:
	od:
od:

return min1:

end:



#GP(L,x,deg) given a set of lists of numbers of the form [[x_1,...,x_n],f(x_1,...,x_n)] and a symbolic variable x returns the polynomial of degree deg in each of the variables which returns f(x_1,..,x_n) when given x_1,...,x_n
GP:=proc(L,x,deg) local i,j,k,num,genpol,Ve2,Ve,T,vars,eqs:

num := nops(L[1][1]):

Ve2 := { seq(convert(i,base,deg+1) ,i=0..(deg+1)^num-1) }:

Ve := {}:

for i from 1 to nops(Ve2) do

	if nops(Ve2[i]) < num then
		Ve := Ve union {[op(Ve2[i]),0$(num-nops(Ve2[i]))]}:
	else
		Ve := Ve union {Ve2[i]}:
	fi:
od:


genpol := 0:

for j from 1 to nops(Ve) do
	genpol := genpol + T(op(Ve[j]))*mul(x[i]^Ve[j][i],i=1..num):
od:

vars := {seq(T(op(Ve[j])),j=1..nops(Ve))}:

eqs:= {seq(subs(seq(x[i]=L[j][1][i],i=1..num),genpol)-L[j][2],j=1..nops(L))}:


return subs(solve(eqs,vars),genpol):


end:

#GR(L,x,deg) given a set of lists of numbers of the form [[x_1,...,x_n],f(x_1,...,x_n)] and a symbolic variable x returns the rational function of degree deg in each of the variables which returns f(x_1,..,x_n) when given x_1,...,x_n
GR:=proc(L,x,deg) local i,j,k,num,genrat,Ve2,Ve,S,T,vars,eqs,genpol1,genpol2:

num := nops(L[1][1]):

Ve2 := { seq(convert(i,base,deg+1) ,i=0..(deg+1)^num-1) }:

Ve := {}:

for i from 1 to nops(Ve2) do

	if nops(Ve2[i]) < num then
		Ve := Ve union {[op(Ve2[i]),0$(num-nops(Ve2[i]))]}:
	else
		Ve := Ve union {Ve2[i]}:
	fi:
od:


genpol1 := 0:
genpol2 := 0:

for j from 1 to nops(Ve) do
	genpol1 := genpol1 + T(op(Ve[j]))*mul(x[i]^Ve[j][i],i=1..num):
	genpol2 := genpol2 + S(op(Ve[j]))*mul(x[i]^Ve[j][i],i=1..num):
od:

genrat := genpol2/genpol1:

vars := {seq(T(op(Ve[j])),j=1..nops(Ve))} union {seq(S(op(Ve[j])),j=1..nops(Ve))}:

eqs:= {seq(subs(seq(x[i]=L[j][1][i],i=1..num),genrat)-L[j][2],j=1..nops(L))}:


return normal(subs(solve(eqs,vars),genrat)):


end:





#Finds if the expression expn is symmetric in b and c
Issym := proc(expn,b,c) local i,j,gu:

gu := expand(expn):

for i from 0 to nops(gu) do
	for j from 0 to i do

		if coeff(coeff(gu,b,i),c,j) <> coeff(coeff(gu,b,j),c,i) then
			print(i,j):
			return 0:
		fi:
	od:
od:

return 1:

end:







##################Functions from GuessHolo2: findrec() and Yafe()
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a,dim:
option remember:

dim := LinearAlgebra[Dimension](f);
if (1+DEGREE)*(1+ORDER)+5+ORDER>dim then
error "Insufficient data for a recurrence of order %1 degree %2",ORDER, DEGREE:
fi:
ope:=0:
var:={}:

ope := add( add(a[i,j]*n^j*N^i, j=0..DEGREE), i=0..ORDER);
var := indets(ope) minus {n,N};

eq := {seq(
  add(subs(n=n0,coeff(ope,N,i))*f[n0+i], i=0..ORDER),
  n0 = 1..dim-ORDER)};

var1:=solve(eq,var):

kv := map(lhs, select( x->evalb(op(1,x) = op(2,x)), var1));

ope:=eval(ope, var1);
if ope=0 then return FAIL end if;

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

`if`( nops(ope)>=1,Yafe(ope[1],N)[2],FAIL);
end:



Yafe:=proc(ope,N) local i,ope1,coe1,L:
    if ope=0 then return (1,0) end if;
    ope1:=expand(ope):
    L:=degree(ope1,N):
    coe1:=coeff(ope1,N,L):
    ope1:=normal(ope1/coe1):
    # ope1:=normal(ope1):
    ope1 := collect(ope1, N, factor);
    factor(coe1),ope1:
end: