print(`This is the package RS`):
print(`written by Arvind Ayyer. Version of 30 November, 2008.`):
print(`It is written for the analysis of FPLs,`):
print(`particularly the Razumov-Stroganov conjecture.`):
print(`FPLs are stored in an internal format in two-dimensional arrays.`):
print(`To see the list of procedures, type Help();`):
print(`For help on a specific procedure, type Help(procedure_name);`):


Help := proc()

if args=NULL then

print(`Contains the following procedures:`):
print(`ASM: nasm,asm,asmw1,asmtofpl,fpltoasm,fpltograph.`):
print(`FPL: fpls,A,printfpl,fplconn,npitofpl,pitofpl,ei,invei.`):
print(`RS Conjecture: revrs,testrs,rsmatrix,redrsmatrix`):
print(`Alternating Paths: findpath,apsq,appath,newaltfpl`):
print(``):
print(`WARNING: The lattice notation is unconventional. The origin is on the top left. The -y axis coordinate is first followed by the +x axis coordinate`):

fi:

if nops([args])=1 and op(1,[args])=`nasm` then
print(`nasm(n) returns the number of ASMs of size n`):
print(`For example, try nasm(3);`):
fi:


if nops([args])=1 and op(1,[args])=`asm` then
print(`asm(n) returns all ASMs of size n`):
print(`For example, try asm(3);`):
fi:


if nops([args])=1 and op(1,[args])=`asmw1` then
print(`asmw1(n,k) returns all ASMs of size n whose 1`):
print(`in the first row is in the k'th column.`):
print(`For example, try asm(3,2);`):
fi:

if nops([args])=1 and op(1,[args])=`fpltograph` then
print(`fpltograph(fpl) for an fpl in the internal formal returns `):
print(`the same fpl as a graph in the subset of the 2d lattice`):
print(`For example, try fpltograph([[5, 3, 2], [3, 2, 5], [2, 5, 3]]);`):
fi:

if nops([args])=1 and op(1,[args])=`asmtofpl` then
print(`asmtofpl(asm) returns the fpl corresponding`):
print(`to the given asm.`):
print(`For example, try asmtofpl([[0, 1, 0], [0, 0, 1], [1, 0, 0]]);`):
fi:


if nops([args])=1 and op(1,[args])=`fpltoasm` then
print(`asmtofpl(fpl) returns the ASM corresponding`):
print(`to the given fpl (in the internal format).`):
print(`For example, try fpltoasm([[5, 1, 4], [3, 6, 1], [2, 5, 3]]);`):
fi:


if nops([args])=1 and op(1,[args])=`fpls` then
print(`fpls(n) returns a list of link patterns of size 2n`):
print(`along with the number of fpls in that link pattern.`):
print(`For example, try fpls(2);`):
fi:

if nops([args])=1 and op(1,[args])=`printfpl` then
print(`printfpl(fpl) prints the fpl in a screen readable format.`):
print(`The 6 possible connections are drawn as ||,==,-|,|-,_|,|_.`):
print(`For example, try printfpl([[5, 1, 4], [3, 6, 1], [2, 5, 3]]);`):
fi:

if nops([args])=1 and op(1,[args])=`fplconn` then
print(`fplconn(fpl) returns the link pattern of that fpl.`):
print(`For example, try fplconn([[5, 1, 4], [3, 6, 1], [2, 5, 3]]);`):
fi:

if nops([args])=1 and op(1,[args])=`npitofpl` then
print(`npitofpl(pi) given a link pattern returns all fpls`):
print(`which have the same link pattern (in the internal format).`):
print(`For example, try npitofpl({{1, 6}, {2, 3}, {4, 5}});`):
fi:


if nops([args])=1 and op(1,[args])=`pitofpl` then
print(`npitofpl(pi) given a link pattern returns all fpls`):
print(`which have the same link pattern in the form of`):
print(`edges in two dimensional square lattice.`):
print(`For example, try pitofpl({{1, 6}, {2, 3}, {4, 5}});`):
fi:


if nops([args])=1 and op(1,[args])=`ei` then
print(`ei(pi,i) given a link pattern of size n and an integer from `):
print(`1 to 2n return the new link pattern got by connecting i to i+1`):
print(`cyclically and their partners to each other.`):
print(`For example, try ei({{1, 6}, {2, 3}, {4, 5}},1);`):
fi:

if nops([args])=1 and op(1,[args])=`invei` then
print(`invei(pi) given a link pattern returns all other link patterns pi'`):
print(`that return pi under the action of ei(pi',i) along with the  i's.`):
print(`For example, try invei({{1, 6}, {2, 3}, {4, 5}});`):
fi:

if nops([args])=1 and op(1,[args])=`A` then
print(`A(pi) given a link patterns returns the number of FPLs that`):
print(`have that particular link pattern.`):
print(`For example, try A({{1, 2}, {3, 4}});`):
fi:

if nops([args])=1 and op(1,[args])=`testrs` then
print(`testrs(n) tests the RS conjecture for fpls of size n.`):
print(`If true, you should get a list of zeros of size Catalan(n)`):
print(`For example, try testrs(2);`):
fi:

if nops([args])=1 and op(1,[args])=`revrs` then
print(`revrs(n) assumes the RS conjecture to find the number of FPLs`):
print(`for each link pattern. Faster that fpls(n).`):
print(`For example, try revrs(2);`):
fi:

if nops([args])=1 and op(1,[args])=`rsmatrix` then
print(`rsmatrix(n) returns the matrix in the RS conjecture in the `):
print(`link pattern basis ordered in no particular way.`):
print(`For example, try rsmatrix(3);`):
fi:

if nops([args])=1 and op(1,[args])=`redrsmatrix` then
print(`redrsmatrix(n) returns the reduced matrix in the RS conjecture in the `):
print(`basis in which link patterns are considered equivalent upto rotation.`):
print(`For example, try redrsmatrix(3);`):
fi:


if nops([args])=1 and op(1,[args])=`findpath` then
print(`findpath(fpl,i) given an fpl and an endpoint i finds the path `):
print(`on the lattice that leads from i to its endpoint.`):
print(`For example, try findpath([[5, 1], [1, 3]],1);`):
fi:

if nops([args])=1 and op(1,[args])=`apsq` then
print(`apsq(fpl,sq) given an fpl and a square whose top left corner is (i,j)`):
print(`returns all  alternating paths in the fpl which touch atleast two `):
print(`of the four possible edges in that square.`):
print(`For example, try apsq([[5, 3, 2], [1, 4, 5], [6, 1, 3]],[1,1]);`):
fi:

if nops([args])=1 and op(1,[args])=`appath` then
print(`appath(fpl,path) given an fpl and a list of vertices forming part of an`):
print(`alternating path, returns all complete alternating paths which contain it.`):
print(`For example, try appath([[5, 3, 2], [1, 4, 5], [6, 1, 3]],[[1,1],[2,1]]);`):
fi:

if nops([args])=1 and op(1,[args])=`newaltfpl` then
print(`newaltfpl(fpl,ap) given an fpl and a complete alternating path returns`):
print(`the new fpl obtained by the involution.`):
print(`For example, try newaltfpl([[5, 1, 4], [1, 1, 1], [6, 1, 3]],[[1, 1], [2, 1], [2, 2], [1, 2], [1, 1]]);`):
fi:

if nops([args])=1 and op(1,[args])=`` then
print(` `):
print(``):
print(`For example, try ;`):
fi:

end:


##################################ASM procedures############################
nasm := proc(n) local k:

return mul((3*k+1)!/(n+k)!,k=0..n-1):

end:

#asm(n) returns the ASMs of size n
asm := proc(n) local gogs, asms:

gogs := GOGset(n,n):

asms := {seq(op(GOGTOASM(op(i,gogs))),i=1..nops(gogs))}:

return asms:

end:

#asmw1(n,k) returns the ASMs with the 1 in the first row in the kth column
asmw1:= proc(n,k) local i,j,gogs,asms,asmk:

if k<0 or k>n then
	ERROR(`The column where the 1 should be located is not correct`):
fi:

gogs := GOGset(n,n):

asms := {seq(op(GOGTOASM(op(i,gogs))),i=1..nops(gogs))}:

asmk := {}:

for i from 1 to nops(asms) do
	if asms[i][1,k] = 1 then
		asmk := asmk union {asms[i]}:
	fi:
od:

return asmk:

end:



printfpl:=proc(cnt) local i,j:

for i from 1 to nops(cnt) do
	printf("\t"):
	for j from 1 to nops(cnt) do
		if cnt[i,j]=1 then
			printf("||\t"):
		elif cnt[i,j]=2 then
			printf("==\t"):
		elif cnt[i,j]=3 then
			printf("|_\t"):
		elif cnt[i,j]=4 then
			printf("|-\t"):
		elif cnt[i,j]=5 then
			printf("-|\t"):
		elif cnt[i,j]=6 then
			printf("_|\t"):
		else
			printf("00\t"):
		fi:
	od:
	printf("\n"):
od:

end:


fztop := proc(cnt,i,j) local k:

if i=1 then
	if j=1 then
		return 5:
	elif j mod 2 = 1 then
		if cnt[i,j-1] = 2 or cnt[i,j-1]=3 or cnt[i,j-1]=4 then
			return 5:
		else
			return 4:
		fi:
	elif j mod 2 = 0 then
		if cnt[i,j-1] = 2 or cnt[i,j-1]=3 or cnt[i,j-1]=4 then
			return 6:
		else
			return 3:
		fi:
	fi:
else
	if cnt[i,j-1] = 2 or cnt[i,j-1] = 3 or cnt[i,j-1] = 4 then
		if cnt[i-1,j] = 1 or cnt[i-1,j] = 4 or cnt[i-1,j] = 5 then	
			return 6:
		else
			return 5:
		fi:
	else
		if cnt[i-1,j] = 1 or cnt[i-1,j] = 4 or cnt[i-1,j] = 5 then	
			return 3:
		else
			return 4:
		fi:
	fi:
fi:


end:

fzright := proc(cnt,i,j) local k:

if j=nops(cnt) then
	if i=1 and nops(cnt) mod 2 = 0 then
		return 6:
	elif i=1 and nops(cnt) mod 2 = 1 then
		return 4:
	elif (i mod 2 = 1 and nops(cnt) mod 2 = 1) or (i mod 2 = 0 and nops(cnt) mod 2 = 0)then
		if cnt[i-1,j] = 2 or cnt[i-1,j]=3 or cnt[i-1,j]=6 then
			return 4:
		else
			return 3:
		fi:
	elif (i mod 2 = 1 and nops(cnt) mod 2 = 0) or (i mod 2 = 0 and nops(cnt) mod 2 = 1)then
		if cnt[i-1,j] = 2 or cnt[i-1,j]=3 or cnt[i-1,j]=6 then
			return 5:
		else
			return 6:
		fi:
	fi:
else
	if cnt[i,j+1] = 2 or cnt[i,j+1] = 5 or cnt[i,j+1] = 6 then
		if cnt[i-1,j] = 1 or cnt[i-1,j] = 4 or cnt[i-1,j] = 5 then	
			return 3:
		else
			return 4:
		fi:
	else
		if cnt[i-1,j] = 1 or cnt[i-1,j] = 4 or cnt[i-1,j] = 5 then	
			return 6:
		else
			return 5:
		fi:
	fi:

fi:


end:

fzbot := proc(cnt,i,j) local k:

if i=nops(cnt) then
	if j=nops(cnt) then
		return 3:
	elif (nops(cnt)-j) mod 2 = 1 then
		if cnt[i,j+1] = 1 or cnt[i,j+1]=3 or cnt[i,j+1]=4 then
			return 5:
		else
			return 4:
		fi:
	elif (nops(cnt)-j) mod 2 = 0 then
		if cnt[i,j+1] = 1 or cnt[i,j+1]=3 or cnt[i,j+1]=4 then
			return 6:
		else
			return 3:
		fi:
	fi:
else
	if cnt[i,j+1] = 2 or cnt[i,j+1] = 5 or cnt[i,j+1] = 6 then
		if cnt[i+1,j] = 1 or cnt[i+1,j] = 3 or cnt[i+1,j] = 6 then	
			return 4:
		else
			return 3:
		fi:
	else
		if cnt[i+1,j] = 1 or cnt[i+1,j] = 3 or cnt[i+1,j] = 6 then	
			return 5:
		else
			return 6:
		fi:
	fi:

fi:


end:


fzleft := proc(cnt,i,j) local k:

if j=1 then
	if i=nops(cnt) and nops(cnt) mod 2 = 1 then
		return 6:
	elif i=nops(cnt) and nops(cnt) mod 2 = 0 then
		return 4:
	elif i mod 2 = 1 then
		if cnt[i+1,j] = 2 or cnt[i+1,j]=4 or cnt[i+1,j]=5 then
			return 6:
		else
			return 5:
		fi:
	elif i mod 2 = 0 then
		if cnt[i+1,j] = 2 or cnt[i+1,j]=4 or cnt[i+1,j]=5 then
			return 3:
		else
			return 4:
		fi:
	fi:
else
	if cnt[i,j-1] = 2 or cnt[i,j-1] = 3 or cnt[i,j-1] = 4 then
		if cnt[i+1,j] = 1 or cnt[i+1,j] = 3 or cnt[i+1,j] = 6 then	
			return 5:
		else
			return 6:
		fi:
	else
		if cnt[i+1,j] = 1 or cnt[i+1,j] = 3 or cnt[i+1,j] = 6 then	
			return 4:
		else
			return 3:
		fi:
	fi:

fi:


end:




#asmtofpl(asm) returns and draws an FPL for a given ASM
asmtofpl := proc(asm) local i,j,k,dim,cnt:

dim := nops(convert(asm,listlist));

#cnt will mark each corner according to the type of vertex it has
#cnt(||)=1, cnt(__)=2, cnt(|_)=3, cnt(|-)=4, cnt(-|)=5, cnt(_|)=6
cnt:=[[0$dim]$dim]:

for i from 1 to ceil(dim/2) do

	for j from i to dim+1-i do
		
		#Top Row
		if cnt[i,j]=0 and (i+j) mod 2 = 1 then
			if asm[i,j]=0 then
				cnt[i,j]:=fztop(cnt,i,j):
			elif asm[i,j]=1 then
				cnt[i,j]:=1:
			else
				cnt[i,j]:=2:
			fi:
		elif cnt[i,j]=0 and (i+j) mod 2 = 0 then
			if asm[i,j]=0 then
				cnt[i,j]:=fztop(cnt,i,j):
			elif asm[i,j]=1 then
				cnt[i,j]:=2:
			else
				cnt[i,j]:=1:
			fi:
		fi:

		#Right column
		if cnt[j,dim+1-i]=0 and (dim+1-i+j) mod 2 = 1 then
			if asm[j,dim+1-i]=0 then
				cnt[j,dim+1-i]:=fzright(cnt,j,dim+1-i):
			elif asm[j,dim+1-i]=1 then
				cnt[j,dim+1-i]:=1:
			else
				cnt[j,dim+1-i]:=2:
			fi:
		elif cnt[j,dim+1-i]=0 and (dim+1-i+j) mod 2 = 0 then
			if asm[j,dim+1-i]=0 then
				cnt[j,dim+1-i]:=fzright(cnt,j,dim+1-i):
			elif asm[j,dim+1-i]=1 then
				cnt[j,dim+1-i]:=2:
			else
				cnt[j,dim+1-i]:=1:
			fi:
		fi:

		#Bottom Row
		if cnt[dim+1-i,dim+1-j]=0 and (2*dim+2-i-j) mod 2 = 1 then
			if asm[dim+1-i,dim+1-j]=0 then
				cnt[dim+1-i,dim+1-j]:=fzbot(cnt,dim+1-i,dim+1-j):
			elif asm[dim+1-i,dim+1-j]=1 then
				cnt[dim+1-i,dim+1-j]:=1:
			else
				cnt[dim+1-i,dim+1-j]:=2:
			fi:
		elif cnt[dim+1-i,dim+1-j]=0 and (2*dim+2-i-j) mod 2 = 0 then
			if asm[dim+1-i,dim+1-j]=0 then
				cnt[dim+1-i,dim+1-j]:=fzbot(cnt,dim+1-i,dim+1-j):
			elif asm[dim+1-i,dim+1-j]=1 then
				cnt[dim+1-i,dim+1-j]:=2:
			else
				cnt[dim+1-i,dim+1-j]:=1:
			fi:
		fi:

		#Left Column
		if cnt[dim+1-j,i]=0 and (i+dim+1-j) mod 2 = 1 then
			if asm[dim+1-j,i]=0 then
				cnt[dim+1-j,i]:=fzleft(cnt,dim+1-j,i):
			elif asm[dim+1-j,i]=1 then
				cnt[dim+1-j,i]:=1:
			else
				cnt[dim+1-j,i]:=2:
			fi:
		elif cnt[dim+1-j,i]=0 and (i+dim+1-j) mod 2 = 0 then
			if asm[dim+1-j,i]=0 then
				cnt[dim+1-j,i]:=fzleft(cnt,dim+1-j,i):
			elif asm[dim+1-j,i]=1 then
				cnt[dim+1-j,i]:=2:
			else
				cnt[dim+1-j,i]:=1:
			fi:
		fi:


	od:
od:

#printcnt(cnt):
return cnt:

end:


#fpltoasm(fpl) returns the asm corresponding to the given fpl
fpltoasm := proc(fpl) local asm,i,j,dim:

dim := nops(fpl):

asm := [[0$dim]$dim]:

for i from 1 to dim do
	for j from 1 to dim do
		if fpl[i,j] > 2 then
			asm[i,j]=0:
		elif fpl[i,j] = 2 then
			if (i+j) mod 2 = 0 then
				asm[i,j]:= 1:
			else
				asm[i,j]:= -1:
			fi:
		elif fpl[i,j] = 1 then
			if (i+j) mod 2 = 0 then
				asm[i,j]:= -1:
			else
				asm[i,j]:= 1:
			fi:
		else
			ERROR(`Something wrong in input`):
		fi:
	od:
od:

return asm:


end:

#asmtogog(asm) given an ASM returns the corresponding monotone triangle
asmtogog := proc(asm) local i,j,k,cnt,gog:

gog := [seq([seq(0,j=1..i)],i=1..nops(asm))]:

for i from 1 to nops(asm) do
	cnt:=1:
	for j from 1 to nops(asm) do
		if add(asm[k,j],k=1..i) = 1 then
			gog[i][cnt]:= j:
			cnt:=cnt+1:
		fi:
	od:
od:

return gog:

end:



##################################FPL procedures############################
#fplconn(fpl) returns the connectivities of the given FPL as a set of sets of two end vertices
fplconn := proc(fpl) local gr,i,j,dim,cc,ends:

dim := nops(fpl):

gr := fpltograph(fpl):

ends := {}:

for i from 1 to dim do
	ends := {seq([2*i-1,0],i=1..ceil(dim/2))} union {seq([0,2*i],i=1..floor(dim/2))}
		union {seq([dim+2-2*i,dim+1],i=1..ceil(dim/2))} union {seq([dim+1,dim+1-2*i],i=1..floor(dim/2))}:
od:

cc:= {seq(CC(gr,ends[i]) intersect ends,i=1..nops(ends))}:

cc :=subs({seq([0,2*i]=i,i=1..floor(dim/2))},cc):
cc :=subs({seq([2*i+dim-2*ceil(dim/2),dim+1]=floor(dim/2)+i,i=1..ceil(dim/2))},cc):
cc :=subs({seq([dim+1,dim+1-2*i]=dim+i,i=1..floor(dim/2))},cc):
cc :=subs({seq([2*ceil(dim/2)+1-2*i,0]=dim+floor(dim/2)+i,i=1..ceil(dim/2))},cc):

return cc:

end:



#fpls(n) returns the number of times a particular link pattern appears for all fpls of size n by n
fpls := proc(n) local i,j,asmn,fpln,lps,uniqlp,numlp:

asmn := asm(n):

fpln := {seq(asmtofpl(asmn[i]),i=1..nops(asmn))}:

lps := [seq(fplconn(fpln[i]),i=1..nops(fpln))]:

uniqlp := {seq(lps[i],i=1..nops(lps))}:

numlp := Array(1..(2*n)!/n!/(n+1)!):

for i from 1 to nops(lps) do
	for j from 1 to nops(uniqlp) do
		if lps[i]=uniqlp[j] then
			numlp[j] := numlp[j]+1:
		fi:
	od:
od:

return {seq([uniqlp[i],numlp[i]],i=1..nops(uniqlp))}:


end:


#npitofpl(pi) calculates all the fpls that contribute to the given pi
npitofpl := proc(pi) local fpln,asmn,lps,i,j,pifpl:

pifpl := []:

asmn := asm(nops(pi)):

fpln := [seq(asmtofpl(asmn[i]),i=1..nops(asmn))]:

lps := [seq(fplconn(fpln[i]),i=1..nops(fpln))]:

for i from 1 to nops(fpln) do
	if pi = lps[i] then
#		printf("%d\n",i):
#		printfpl(fpln[i]):
		pifpl := [op(pifpl),fpln[i]]:
	fi:
od:

return pifpl:

end:




#A(pi) returns the number of FPLs with the given endpoint connectivities pi
A := proc(pi) local i,num,fpla:

fpla := revrs(nops(pi)):

for i from 1 to nops(fpla) do
	if fpla[i][1] = pi then
		num := fpla[i][2]:
	fi:
od:

return num:

end:


#ei(pi,i) returns the new connectivity obtained by the action of e_i(pi)
ei := proc(pi,i) local j,k,npi,ipos,ipluspos,l1,l2,iplus1:

if i<2*nops(pi) then
	iplus1 := i+1:
else
	iplus1 := 1:
fi:

for j from 1 to nops(pi) do
	if pi[j][1] = i then
		ipos := j:
		l1 := pi[j][2]:
	fi:

	if pi[j][2] = i then
		ipos := j:
		l1 := pi[j][1]:
	fi:

	if pi[j][1] = iplus1 then
		ipluspos := j:
		l2 := pi[j][2]:
	fi:

	if pi[j][2] = iplus1 then
		ipluspos := j:
		l2 := pi[j][1]:
	fi:
od:

if ipos < ipluspos then
	npi := {seq(pi[j],j=1..ipos-1),{i,iplus1},seq(pi[j],j=ipos+1..ipluspos-1),{l1,l2},seq(pi[j],j=ipluspos+1..nops(pi))}:
elif ipluspos < ipos then
	npi := {seq(pi[j],j=1..ipluspos-1),{i,iplus1},seq(pi[j],j=ipluspos+1..ipos-1),{l1,l2},seq(pi[j],j=ipos+1..nops(pi))}:
else
	npi := pi:
fi:

return npi:

end:


#invei(pi) given a connectivity pi, returns the set of i's and sigma's such that ei(sigma,i)=pi
invei := proc(pi) local i,j,nc,gu:

nc := numchords(1,2*nops(pi)):

gu:= {}:

for i from 1 to nops(nc) do
	for j from 1 to 2*nops(pi) do
		if ei(nc[i],j) = pi then
			gu := gu union {[j,nc[i]]}:
		fi:
	od:
od:

return gu:

end:

###########################RS Conjecture procedures#############################
#testrs(n) tests the RS conjecture for all configurations in LP(n)
testrs := proc(n) local i,j,k,fpla,S,api,numpi:

fpla := fpls(n):

S := []:
for i from 1 to nops(fpla) do

	api := fpla[i][1]:
	numpi := 2*n*fpla[i][2]:
	for j from 1 to nops(fpla) do
		for k from 1 to 2*n do

			if ei(fpla[j][1],k) = api then
				numpi := numpi - fpla[j][2]:
			fi:
		od:
	od:

	S := [op(S),numpi]:

od:

return S:


end:


Eq1:=proc(pi,n,A)  local i,sig,lu,gu:
gu:=2*n*A[pi]:
lu:=numchords(1,2*n):
for sig in lu do
for i from 1 to 2*n do
 if ei(sig,i)=pi then
  gu:=gu-A[sig]:
 fi:
od:
od:
gu:
end:


#revrs(n) reverses the RS conjecture to solve for A(pi) for all pi's.
revrs := proc(n) local gu,pi,eq,var,pi0,A,i,mu,m:

pi0:={seq({2*i,2*i-1},i=1..n)}:
gu:=numchords(1,2*n):
eq:= {seq(Eq1(pi,n,A), pi in gu)}:
#return eq:

eq:=subs({A[pi0]=1},eq):

var:= {seq(A[pi], pi in gu minus {pi0})} :
var:=solve(eq,var):

mu:={[pi0,1],seq([pi,subs(var,A[pi])],pi in gu minus {pi0} )}:

m:=min(seq(mu[i][2],i=1..nops(mu))):


{seq([mu[i][1],mu[i][2]/m],i=1..nops(mu))}:


end:

#rsmatrix(n) finds the matrix involved in the RS conjecture
rsmatrix := proc(n) local i,j,M,C,A,eqs:

C := numchords(1,2*n):

eqs := {seq(Eq1(C[i],n,A),i=1..nops(C))}:
M := [[0$nops(C)]$nops(C)]:

for i from 1 to nops(eqs) do
	for j from 1 to nops(C) do
		M[i][j] := coeff(eqs[i],A[C[j]]):
	od:
od:

return convert(M,matrix):

end:


#redrsmatrix(n) finds the matrix involved in the RS conjecture when only inequivalent patters are considered
redrsmatrix := proc(n) local i,j,k,l,M,C1,C2,A,eqs,chart:

C1 := numchords(1,2*n):
C2 := ineqchords(n):

chart := {}:

for i from 1 to nops(C1) do
	for j from 1 to nops(C2) do
		for l from 1 to 2*n-1 do
			if subs({seq(k=k+l,k=1..2*n-l),seq(k=k+l-2*n,k=2*n-l+1..2*n)},C1[i]) =C2[j] then
				chart := chart union {[C1[i],C2[j]]}:
			fi:
		od:
	od:
od:

#return chart:


eqs := {seq(Eq1(C2[i],n,A),i=1..nops(C2))}:
eqs := subs({seq(chart[i][1]=chart[i][2],i=1..nops(chart))},eqs):

#return eqs:
M := Array(1..nops(C2),1..nops(C2)):
#M := [[0$nops(C2)]$nops(C2)]:

for i from 1 to nops(eqs) do
	for j from 1 to nops(C2) do
		M[i,j] := coeff(eqs[i],A[C2[j]]):
	od:
od:

return convert(M,matrix):

end:





######################################Procedures about the G_S operator in Wieland's paper#########



#gsonfpl(fpl,i,j) applies the operator G_S defined in wieland's paper on the square with lower left vertex=(i,j)
gsonfpl := proc(fpl,i,j) local n,nfpl,f1,f2,f3,f4:

n := nops(fpl):

if (i<0 or i>n) or (j<0 or j>n) then
	ERROR(`i and j must be between 0 and`,n):
fi:

if i=0 or i=n or j=0 or j=n then
	return fpl:
fi:

nfpl := fpl:

f1 := nfpl[i][j]:
f3 := nfpl[i+1][j]:
f2 := nfpl[i][j+1]:
f4 := nfpl[i+1][j+1]:

#Horizontal lines -> Vertical lines
if (f1=3 or f1=2) and (f2=2 or f2=6) and (f3=2 or f3=4) and (f4=2 or f4=5) then
	if f1=2 then
		nfpl[i][j] := 5:
	else
		nfpl[i][j] := 1:
	fi:

	if f2=2 then
		nfpl[i][j+1] := 4:
	else
		nfpl[i][j+1] := 1:
	fi:

	if f3=2 then
		nfpl[i+1][j] := 6:
	else
		nfpl[i+1][j] := 1:
	fi:

	if f4=2 then
		nfpl[i+1][j+1] := 3:
	else
		nfpl[i+1][j+1] := 1:
	fi:
fi:

# Vertical lines -> Horizontal lines
if (f1=1 or f1=5) and (f2=1 or f2=4) and (f3=1 or f3=6) and (f4=1 or f4=3) then
	if f1=1 then
		nfpl[i][j] := 3:
	else
		nfpl[i][j] := 2:
	fi:

	if f2=1 then
		nfpl[i][j+1] := 6:
	else
		nfpl[i][j+1] := 2:
	fi:

	if f3=1 then
		nfpl[i+1][j] := 4:
	else
		nfpl[i+1][j] := 2:
	fi:

	if f4=1 then
		nfpl[i+1][j+1] := 5:
	else
		nfpl[i+1][j+1] := 2:
	fi:
fi:

#printfpl(nfpl):
return nfpl:


end:

#gsodd(fpl) acts the gs operator on the fpl for all odd squares
gsodd := proc(fpl) local i,j,nfpl,n:

n := nops(fpl):
nfpl := fpl:

for i from 0 to n do
	for j from 0 to n do
		if (i+j) mod 2 = 1 then
			nfpl := gsonfpl(nfpl,i,j):
		fi:
	od:
od:

return nfpl:

end:

#gseven(fpl) acts the gs operator on the fpl for all even squares
gseven := proc(fpl) local i,j,nfpl,n:

n := nops(fpl):
nfpl := fpl:

for i from 0 to n do
	for j from 0 to n do
		if (i+j) mod 2 = 0 then
			nfpl := gsonfpl(nfpl,i,j):
		fi:
	od:
od:

return nfpl:

end:

#gsgraph(n) generates the graph whose vertices are all n by n with directed edges marked (i,j) from f to g if the G_S operator on square (i,j) on f takes it to g
gsgraph := proc(n) local i,j,k,asmn,V,E:

asmn := asm(n):

V := {seq(asmtofpl(asmn[i]),i=1..nops(asmn))}:

E := {}:

for i from 1 to nops(V) do
	for j from 1 to n-1 do
		for k from 1 to n-1 do
			if gsonfpl(V[i],j,k) <> V[i] then
				E := E union {[{V[i],gsonfpl(V[i],j,k)},[j,k]]}:
			fi:
		od:
	od:
od:


return V,E:

end:


#gsmoves(fpl,k) returns all the fpls which can be reached by k moves in the corresponding gsgraph()
gsmoves := proc(fpl,k) local n,S,T,i,j,E:
option remember:

if k=0 then
	return {[fpl,[]]}:
fi:

n := nops(fpl):
E := gsgraph(n)[2]:

S := gsmoves(fpl,k-1):
T := {}:
for i from 1 to nops(S) do
	for j from 1 to nops(E) do
		if member(S[i][1],E[j][1]) then
			T := T union {[op(E[j][1] minus {S[i][1]}),[op(S[i][2]),E[j][2]]]}:
		fi:
	od:
od:

S := T:
T := {}:
for i from 1 to nops(S) do
	if nops(S[i][2]) = nops(convert(S[i][2],set)) then
		T := T union {S[i]}:
	fi:
od:

return T:

end:
	
#printgsmoves(fpl,k) prints the fpls that arise from fpl by k G_S moves
printgsmoves := proc(fpl,k) local s,S:

S := gsmoves(fpl,k):

print(printfpl(fpl),` gives rise to`):
for s in S do
	print(printfpl(s[1]),` by moves`, s[2]):
od:

end:

#gsfpltofpl(f1,f2,n) returns the minimal set of steps needed to get from f1 to f2
gsfpltofpl := proc(f1,f2,n) local i,j,g1,g2:

#n := nops(f1):

for i from 0 to n do
	g1 := [op(gsmoves(f1,i))]:
	g2 := [seq(g1[i][1],i=1..nops(g1))]:
	for j from 1 to nops(g2) do
		if g2[j]=f2 then
			return g1[j][2]:
		fi:
	od:
od:

return FAIL:

end:



#findbij(conn) given a connectivity, lists its fpls, lists the fpls that it is supposed to go to and finds all the square moves that takes the former to the latter
findbij := proc(conn) local i,j,fpl1,fpl2,c2,S:

fpl1 := npitofpl(conn):

c2 := invei(conn):
c2 := {seq(c2[i][2],i=1..nops(c2))}:

fpl2 := {seq(op(npitofpl(c2[i])),i=1..nops(c2))}:


#return fpl2:

S := {}:

for i from 1 to nops(fpl1) do
	for j from 1 to nops(fpl2) do
#		print([fpl1[i],fpl2[j],gsfpltofpl(fpl1[i],fpl2[j])]):
		S := S union {[fpl1[i],fpl2[j],gsfpltofpl(fpl1[i],fpl2[j],3)]}:
	od:
od:

print(fpl1):
#print(fpl2):

return S:

end:

######################################Procedures about Alt paths in FPLs#########

#fpltograph(fpl) returns the special multi-component graph represented by the fpl
fpltograph := proc(fpl) local ve,ed,i,j,dim:

dim := nops(fpl):

ve := {seq(seq([i,j],i=0..dim+1),j=0..dim+1)}:

ed := {}:

for i from 1 to dim do
	for j from 1 to dim do
		if fpl[i,j]=1 then
			ed := ed union {{[i-1,j],[i,j]},{[i,j],[i+1,j]}}:
		elif fpl[i,j]=2 then
			ed := ed union {{[i,j-1],[i,j]},{[i,j],[i,j+1]}}:
		elif fpl[i,j]=3 then
			ed := ed union {{[i-1,j],[i,j]},{[i,j],[i,j+1]}}:
		elif fpl[i,j]=4 then
			ed := ed union {{[i+1,j],[i,j]},{[i,j],[i,j+1]}}:
		elif fpl[i,j]=5 then
			ed := ed union {{[i,j-1],[i,j]},{[i,j],[i+1,j]}}:
		elif fpl[i,j]=6 then
			ed := ed union {{[i,j-1],[i,j]},{[i,j],[i-1,j]}}:
		fi:
	od:
od:


return [ve,ed]:

end:

#graphtofpl(gr) returns the fpl represented by the multi-component graph
graphtofpl := proc(gr) local fpl,i,j,k,n,ed,pt,cpt:

n := sqrt(nops(gr[1]))-2:
ed := gr[2]:

fpl := [[0$n]$n]:

for i from 1 to n do
	for j from 1 to n do
		pt := [i,j]:
		cpt := ctpt(ed,pt):
		if evalb([i+1,j] in cpt) and evalb([i-1,j] in cpt) then
			fpl[i][j] := 1:
		elif evalb([i+1,j] in cpt) and evalb([i,j+1] in cpt) then
			fpl[i][j] := 4:
		elif evalb([i+1,j] in cpt) and evalb([i,j-1] in cpt) then
			fpl[i][j] := 5:
		elif evalb([i-1,j] in cpt) and evalb([i,j+1] in cpt) then
			fpl[i][j] := 3:
		elif evalb([i-1,j] in cpt) and evalb([i,j-1] in cpt) then
			fpl[i][j] := 6:
		else
			fpl[i][j] := 2:
		fi:
	od:
od:	

return fpl:

end:

#findpath(fpl,i) given an fpl and a starting point finds the path in the fplgraph starting at point i
findpath := proc(fpl,i) local n,j,k,ed,cnt,ipt,S,ip1:

n := nops(fpl):
ed := fpltograph(fpl)[2]:

if i <= floor(n/2) then
	ipt := [0,2*i]:
	ip1 := [1,2*i]:
elif i <= n then
	ipt := [2*i-n,n+1]:
	ip1 := [2*i-n,n]:
elif i <= n + floor(n/2) then
	ipt := [n+1,3*n+1-2*i]:
	ip1 := [n,3*n+1-2*i]:
else
	ipt := [4*n+1-2*i,0]:
	ip1 := [4*n+1-2*i,1]:
fi:

S := [ipt,ip1]:
cnt := 1:

for j from 1 while cnt = 1 do
	cnt := 0:
	for k from 1 to nops(ed) do
		if evalb(S[nops(S)] in ed[k]) and not evalb(S[nops(S)-1] in ed[k]) then
			S := [op(S),op(ed[k] minus {S[nops(S)]})]:
			cnt := 1:
		fi:
	od:
od:

return S:

end:

#distpaths(p1,p2) given two paths in an fpl finds the shortest distance points between them
distpaths := proc(p1,p2) local i,j,d,S:

d := 2:
S := []:
for i from 1 to nops(p1) do
	for j from 1 to nops(p2) do
		if abs(p1[i][1]-p2[j][1])+abs(p1[i][2]-p2[j][2]) < d then
			S := [[p1[i],p2[j]]]:
			d := abs(p1[i][1]-p2[j][1])+abs(p1[i][2]-p2[j][2]):
		elif abs(p1[i][1]-p2[j][1])+abs(p1[i][2]-p2[j][2]) = d then
			S := [op(S),[p1[i],p2[j]]]:
		fi:
	od:
od:

return S:

end:

#distpaths4(p1,p2) given two paths in an fpl finds the shortest distance points between them
distpaths4 := proc(p1,p2) local i,j,d,S,f:

d := 4:
S := []:
for i from 2 to nops(p1)-2 do
	for j from 2 to nops(p2)-2 do
		f := abs(p1[i][1]-p2[j][1])+abs(p1[i][2]-p2[j][2])
		  + abs(p1[i+1][1]-p2[j+1][1])+abs(p1[i+1][2]-p2[j+1][2]):
		if f <= d then
			S := [op(S),[[p1[i],p1[i+1]],[p2[j],p2[j+1]]]]:
		fi:
	od:
od:

return S:

end:

#distpaths2(p1,p2) given two paths in an fpl finds the shortest distance points between them
distpaths2 := proc(p1,p2) local i,j,d,S,f:

d := 4:
S := []:
for i from 2 to nops(p1)-2 do
	for j from 2 to nops(p2)-2 do
		f := abs(p1[i][1]-p2[j][1])+abs(p1[i][2]-p2[j][2])
		  + abs(p1[i+1][1]-p2[j+1][1])+abs(p1[i+1][2]-p2[j+1][2]):
		if f < d then
			S := [[[p1[i],p1[i+1]],[p2[j],p2[j+1]]]]:
			d := f:
		elif f = d then
			S := [op(S),[[p1[i],p1[i+1]],[p2[j],p2[j+1]]]]:
		fi:
	od:
od:

return S:

end:


#ctpt(ed,pt) given the edges in the fplgraph and a lattice point, finds the other two points connected to it in the fpl
ctpt := proc(ed,pt) local a,b,S:

S := {}:

for a from 1 to nops(ed) do
	if evalb(pt in ed[a]) then
		S := S union {op(ed[a] minus {pt})}:
	fi:
od:

return S:

end: 

#looparea(L) given a loop in the integer lattice finds the area enclosed by the loop
looparea := proc(L) local i,j,minx,maxx,miny,maxy:
option remember:

if nops(L)=0 then
	return 0:
fi:

minx := min(seq(L[i][1],i=1..nops(L))):
maxx := max(seq(L[i][1],i=1..nops(L))):
miny := min(seq(L[i][2],i=1..nops(L))):
maxy := max(seq(L[i][2],i=1..nops(L))):

#return minx,miny,maxx,maxy:
for i from 1 to nops(L)-1 do
	for j from i+1 to nops(L)-1 do
		if L[i]=L[j] then
			return looparea([op(1..i,L),op(j+1..nops(L),L)]) + looparea([op(i..j,L)]):
		fi:
	od:
od:


if {seq(op([[maxx,i],[minx,i]]),i=miny..maxy),seq(op([[i,maxy],[i,miny]]),i=minx..maxx)}={op(L)} then
	return (maxx-minx)*(maxy-miny):
fi:


for i from 1 to nops(L)-2 do
	if L[i]=L[i+2] then
		return looparea([op(1..i,L),op(i+3..nops(L),L)]):
	fi:
od:

for i from 1 to nops(L)-2 do
	if L[i][1]=maxx-1 and L[i+1][1]=maxx then
		if L[i+2][2]>L[i+1][2] then
			return 1 + looparea([op(1..i,L),[L[i][1],L[i][2]+1],op(i+2..nops(L),L)]):
		elif L[i+2][2]<L[i+1][2] then
			return 1 + looparea([op(1..i,L),[L[i][1],L[i][2]-1],op(i+2..nops(L),L)]):
		fi:
	fi:
od:

end:

#nloops(fpl) calculates the number of loops in the given fpl
nloops := proc(fpl) local i,j,n,gr,V,path:

n := nops(fpl):
gr := fpltograph(fpl)[2]:

for i from 1 to 2*n do
	path := findpath(fpl,i):
	path := {seq({path[i],path[i+1]},i=1..nops(path)-1)}:
	gr := gr minus path:
od:

if nops(gr)=0 then
	return 0:
fi:

#return gr:
V := {seq(seq(gr[i][j],j=1..2),i=1..nops(gr))}:

return `if`(CCD([V,gr])=0,0,nops(CCD([V,gr]))):

end:


#apsq(fpl,sq) finds all alt paths involving sq in the fpl
apsq := proc(fpl,sq) local i,j,k,b1,w1,S,spt,npt,S1,Sprime,B,W,n:
option remember:

n := nops(fpl):

S := {seq(seq({[j,k],[j+1,k]},j=0..n),k=1..n)}
union {seq(seq({[j,k],[j,k+1]},k=0..n),j=1..n)}:

B := fpltograph(fpl)[2]:

W := S minus B:

#return B,W:

b1 := {}:
w1 := {}:

S := {}:

if evalb({sq,[sq[1]+1,sq[2]]} in B) then
	b1 := b1 union {{sq,[sq[1]+1,sq[2]]}}:
else
	w1 := w1 union {{sq,[sq[1]+1,sq[2]]}}:
fi:

if evalb({sq,[sq[1],sq[2]+1]} in B) then
	b1 := b1 union {{sq,[sq[1],sq[2]+1]}}:
else
	w1 := w1 union {{sq,[sq[1],sq[2]+1]}}:
fi:

if evalb({[sq[1]+1,sq[2]],[sq[1]+1,sq[2]+1]} in B) then
	b1 := b1 union {{[sq[1]+1,sq[2]],[sq[1]+1,sq[2]+1]}}:
else
	w1 := w1 union {{[sq[1]+1,sq[2]],[sq[1]+1,sq[2]+1]}}:
fi:

if evalb({[sq[1],sq[2]+1],[sq[1]+1,sq[2]+1]} in B) then
	b1 := b1 union {{[sq[1],sq[2]+1],[sq[1]+1,sq[2]+1]}}:
else
	w1 := w1 union {{[sq[1],sq[2]+1],[sq[1]+1,sq[2]+1]}}:
fi:

for i from 1 to nops(b1) do
	for j from 1 to nops(w1) do
		if nops(b1[i] intersect w1[j]) > 0 then
			spt := op(b1[i] minus w1[j]):
			npt := op(b1[i] intersect w1[j]):
			S := S union aps1(B,W,[spt,npt,op(w1[j] minus {npt})],spt):
			#S := S union {[spt,npt,op(w1[j] minus {npt})]}:
		fi:
	od:
od:

S1 := S:
S := {}:

for i from 1 to nops(S1) do
	Sprime := {seq({op(S[j])},j=1..nops(S))}:
	if not evalb({op(S1[i])} in Sprime) then
		S := S union {S1[i]}:
	fi:
od:

return S:

end:

#appath(fpl,P) given an fpl and a part of an alternating path on it returns all allowed alternate fpls from it
appath := proc(fpl,P) local n,i,j,S,B,W:

n := nops(fpl):

S := {seq(seq({[j,k],[j+1,k]},j=0..n),k=1..n)}
union {seq(seq({[j,k],[j,k+1]},k=0..n),j=1..n)}:

B := fpltograph(fpl)[2]:
W := S minus B:

S := aps1(B,W,P,P[1]):

return S:
#return seq(newaltfpl(fpl,S[i]),i=1..nops(S)):

end:

#appath1(fpl,P) given an fpl and a part of an alternating path on it returns all allowed alternate fpls from it
appath1 := proc(fpl,P) local n,i,j,S,B,W:

n := nops(fpl):

S := {seq(seq({[j,k],[j+1,k]},j=0..n),k=1..n)}
union {seq(seq({[j,k],[j,k+1]},k=0..n),j=1..n)}:

B := fpltograph(fpl)[2]:
W := S minus B:

S := aps1(B,W,P,P[1]):

return seq([newaltfpl(fpl,S[i]),S[i]],i=1..nops(S)):

end:


#aps1(B,W,stpart,ept) returns all alt paths with starting edge seq stpart and last 
aps1 := proc(B,W,stpart,ept) local i,j,spt,Wpt,Bpt,p1,S:
option remember:

if stpart[nops(stpart)]=ept then
	return {stpart}:
fi:

spt := stpart[nops(stpart)]:
#p1 := {[spt[1]+1,spt[2]],[spt[1]-1,spt[2]],[spt[1],spt[2]+1],[spt[1],spt[2]-1]} minus {stpart[nops(stpart)-1]}:
p1 := {[spt[1]+1,spt[2]],[spt[1]-1,spt[2]],[spt[1],spt[2]+1],[spt[1],spt[2]-1]} minus {op(2..nops(stpart)-1,stpart)}:
S := {}:

#return p1:

if evalb({stpart[nops(stpart)-1],stpart[nops(stpart)]} in B) then
	Wpt := {}:
	for i from 1 to nops(p1) do
		if evalb({spt,p1[i]} in W) then
			Wpt := Wpt union {p1[i]}:
		fi:
	od:
	if nops(Wpt) = 0 then
		return {}:
	else
		for j from 1 to nops(Wpt) do
			S := S union aps1(B,W,[op(stpart),Wpt[j]],ept):
		od:
	fi:
elif evalb({stpart[nops(stpart)-1],stpart[nops(stpart)]} in W) then
	Bpt := {}:
	for i from 1 to nops(p1) do
		if evalb({spt,p1[i]} in B) then
			#print(p1[i]):
			Bpt := Bpt union {p1[i]}:
		fi:
	od:
	if nops(Bpt) = 0 then
		return {}:
	else
		for j from 1 to nops(Bpt) do
			S := S union aps1(B,W,[op(stpart),Bpt[j]],ept):
		od:
	fi:
fi:

end:


#newaltfpl(fpl,AP) given an fpl and an alt path AP finds the new FPLs given by an involution on the alt path
newaltfpl := proc(fpl,AP) local B,n,gr,B2,i,j:

n := nops(fpl):

gr := fpltograph(fpl):
B := gr[2]:

B2 := B:

for i from 1 to nops(AP)-2 do
	if evalb({AP[i],AP[i+1]} in B) then
		B2 := B2 minus {{AP[i],AP[i+1]}} union {{AP[i+1],AP[i+2]}}:
	fi:
od:

#return B2:
return graphtofpl([gr[1],B2]):

end:

#allowedaltfpl(fpl,sq) given an fpl and a square finds all the alt path fpls which involve the square and have the correct connectivity
allowedaltfpl := proc(fpl,sq) local conn,alconn,i,j,altpaths,S,S2,m,c1:

conn := fplconn(fpl):

alconn := {seq(ei(conn,i),i=1..2*nops(fpl))} minus {conn}:
#alconn := invei(conn)
#alconn := {seq(alconn[i][2],i=1..nops(alconn))} minus {conn}:

altpaths := apsq(fpl,sq):

S := {}:

for i from 1 to nops(altpaths) do
	c1 := fplconn(newaltfpl(fpl,altpaths[i])):
	if evalb(c1 in alconn) then
		#S := S union {[altpaths[i],newaltfpl(fpl,altpaths[i])]}:
		S := S union {newaltfpl(fpl,altpaths[i])}:
	fi:
od:

return S:

end:

#listaltfpls(n) given the size of the fpls lists all allowed fpls got by the alt path method with the number of times they appeared
listaltfpls := proc(n) local i,j,fpls,asms,S:

asms := asm(n):
fpls := [seq(asmtofpl(asms[i]),i=1..nops(asms))]:

#S := []:
#for i from 1 to nops(fpls) do
#	for j from 1 to 2*n do
#		if nops(fplalt(fpls[i],j)) = 1 then
#			S := [op(S),op(fplalt(fpls[i],j))]:
#		fi:
#	od:
#od:

S := [seq(seq(fplalt22(fpls[i],j),j=1..2*n),i=1..nops(fpls))]:
return convert(S, multiset):

end:

#listaltfpls3(n) given the size of the fpls lists all allowed fpls got by the alt path method with the number of times they appeared
listaltfpls3 := proc(n) local i,j,fpls,asms,S,S2:

asms := asm(n):
fpls := [seq(asmtofpl(asms[i]),i=1..nops(asms))]:

S := [seq(seq(fplalt22(fpls[i],j),j=1..2*n),i=1..nops(fpls))]:
S2 := [seq([fpls[i],[]],i=1..nops(fpls))]:

for i from 1 to nops(S) do
	for j from 1 to nops(fpls) do
		if S[i][1] = fpls[j] then
			S2[j] := [S2[j][1],[op(S2[j][2]),S[i][2]]]:
		fi:
	od:
od:

return S2:

end:


#detaltfpl(n,fpl) given the size of the fpls and a given fpl detects how it appeared in listaltfpls(n)
detaltfpl := proc(n,fpl) local i,j,fpls,asms,S:

asms := asm(n):
fpls := [seq(asmtofpl(asms[i]),i=1..nops(asms))]:

S := {}:
for i from 1 to nops(fpls) do
	for j from 1 to 2*n do
		if evalb(fpl in fplalt2(fpls[i],j)) then
			S := S union {{fpls[i],j}}:
		fi:
	od:
od:

return S:

end:

#nfplalt(fpl,i) lists naively the alt paths got from starting point i in the fpl
nfplalt := proc(fpl,i) local j,k,S,B,W,n,conn,iconn,ip1,ipt,ipt1,impt,F,P,F2,P2:

n := nops(fpl):
conn := fplconn(fpl):
iconn := ei(conn,i):

if iconn=conn then
	return {fpl}:
fi:

S := {seq(seq({[j,k],[j+1,k]},j=0..n),k=1..n)}
union {seq(seq({[j,k],[j,k+1]},k=0..n),j=1..n)}:

B := fpltograph(fpl)[2]:
W := S minus B:

if i < 2*n then
	ip1:=i+1:
else
	ip1:=1:
fi:

if i <= floor(n/2) then
	ipt := [1,2*i]:
elif i <= n then
	ipt := [2*i-n,n]:
elif i <= n + floor(n/2) then
	ipt := [n,3*n+1-2*i]:
else
	ipt := [4*n+1-2*i,1]:
fi:

if ip1 <= floor(n/2) then
	ipt1 := [1,2*ip1]:
elif ip1 <= n then
	ipt1 := [2*ip1-n,n]:
elif ip1 <= n + floor(n/2) then
	ipt1 := [n,3*n+1-2*ip1]:
else
	ipt1 := [4*n+1-2*ip1,1]:
fi:


if evalb(type((ipt[1]+ipt1[1])/2,integer)) and evalb(type((ipt[2]+ipt1[2])/2,integer)) then
	impt := [(ipt[1]+ipt1[1])/2,(ipt[2]+ipt1[2])/2]:
	if evalb({ipt,impt} in B) then
		F := [appath(fpl,[ipt,impt,ipt1])]:
		P := [op(aps1(B,W,[ipt,impt,ipt1],ipt))]:
	else
		F := [appath(fpl,[ipt1,impt,ipt])]:
		P := [op(aps1(B,W,[ipt1,impt,ipt],ipt1))]:
	fi:
elif (ipt[1]=1 and ipt[2]=n) or (ipt1[1]=1 and ipt1[2]=n) then
	impt := [1,n]:
	if evalb({[1,n],[2,n]} in B) then
		F := [appath(fpl,[[2,n],[1,n],[1,n-1]])]:
		P := [op(aps1(B,W,[[2,n],[1,n],[1,n-1]],[2,n]))]:
	else
		F := [appath(fpl,[[1,n-1],[1,n],[2,n]])]:
		P := [op(aps1(B,W,[[1,n-1],[1,n],[2,n]],[1,n-1]))]:
	fi:
elif (ipt[1]=n and ipt[2]=n) or (ipt1[1]=n and ipt1[2]=n) then
	impt := [n,n]:
	if evalb({[n,n],[n-1,n]} in B) then
		F := [appath(fpl,[[n-1,n],[n,n],[n,n-1]])]:
		P := [op(aps1(B,W,[[n-1,n],[n,n],[n,n-1]],[n-1,n]))]:
	else
		F := [appath(fpl,[[n,n-1],[n,n],[n-1,n]])]:
		P := [op(aps1(B,W,[[n,n-1],[n,n],[n-1,n]],[n,n-1]))]:
	fi:
elif (ipt[1]=1 and ipt[2]=1) or (ipt1[1]=1 and ipt1[2]=1) then
	impt := [1,1]:
	if evalb({[1,1],[2,1]} in B) then
		F := [appath(fpl,[[2,1],[1,1],[1,2]])]:
		P := [op(aps1(B,W,[[2,1],[1,1],[1,2]],[2,1]))]:
	else
		F := [appath(fpl,[[1,2],[1,1],[2,1]])]:
		P := [op(aps1(B,W,[[1,2],[1,1],[2,1]],[1,2]))]:
	fi:
else
	impt := [n,1]:
	if evalb({[n,1],[n,2]} in B) then
		F := [appath(fpl,[[n,2],[n,1],[n-1,1]])]:
		P := [op(aps1(B,W,[[n,2],[n,1],[n-1,1]],[n,2]))]:
	else
		F := [appath(fpl,[[n-1,1],[n,1],[n,2]])]:
		P := [op(aps1(B,W,[[n-1,1],[n,1],[n,2]],[n-1,1]))]:
	fi:
fi:

return F,P;

F2 := []:
P2 := []:
for j from 1 to nops(F) do
	if fplconn(F[j]) = iconn then
		F2 := [op(F2),F[j]]:
		P2 := [op(P2),P[j]]:
	fi:
od:

if nops(F2) = 0 then
	print(fpl,i,`None`):
fi:

return F2:

end:



#fplalt(fpl,i) lists the alt paths got from starting point i in the fpl
fplalt := proc(fpl,i) local j,k,B,W,W2,W3,S3,S1,S2,S,T,T2,n,ipt,ip1,conn,iconn,cnt,d:

n := nops(fpl):
conn := fplconn(fpl):
iconn := ei(conn,i):

if iconn=conn then
	return {fpl}:
fi:

S := {seq(seq({[j,k],[j+1,k]},j=0..n),k=1..n)}
union {seq(seq({[j,k],[j,k+1]},k=0..n),j=1..n)}:

B := fpltograph(fpl)[2]:
W := S minus B:


if i < 2*n then
	ip1:=i+1:
else
	ip1:=1:
fi:

S := distpaths(findpath(fpl,i),findpath(fpl,ip1)):

T := []:
W2 := []:
cnt := 0:
for j from 1 to nops(S) while cnt=0 do
	S1 := {}:
	for k from 1 to nops(B) do
		if evalb(S[j][1] in B[k]) then
			S1 := S1 union B[k] minus {S[j][1]}:
		fi:
	od:
	S2 := [seq(appath(fpl,[S1[k],S[j][1],S[j][2]]),k=1..nops(S1))]:
	S3 := [seq(op(aps1(B,W,[S1[k],S[j][1],S[j][2]],S1[k])),k=1..nops(S1))]:
	for k from 1 to nops(S2) do
		if fplconn(S2[k]) = iconn then
			T := [op(T),S2[k]]:
			W2 := [op(W2),S3[k]]:
			cnt := 1:
		fi:
	od:
	if nops(T) > 1 then
		d := min(seq(nops(W2[k]),k=1..nops(W2))):
		W3 := []:
		T2 := []:
		for k from 1 to nops(W2) do
			if nops(W2[k]) = d then
				W3 := [op(W3),W2[k]]:
				T2 := [op(T2),T[k]]:
			fi:
		od:
		T := T2:
		if nops(W3)>1 then
			print(fpl,i,W3):
		fi:
	fi:
od:

return T:

end:


inv2 := proc(fpl) local i,j,asms,fpl1,S,wfpl,n:

wfpl := {op(wmapfpl(fpl))}:
n := nops(fpl):
asms := asm(n):
fpl1 := [seq(asmtofpl(asms[i]),i=1..nops(asms))]:
S := []:
for i from 1 to nops(fpl1) do
	if evalb(op(fplalt2(fpl1[i],2*n))[1] in wfpl) then
		S := [op(S),fpl1[i]]:
	fi:
od:

return S:

end:

#fplalt22(fpl,i) uses wmapfpl() to determine the alt fpl
fplalt22 := proc(fpl,i) local n,j,k,wfpl,nfpl:

n := nops(fpl):
if i=2*n then
	return fplalt3(fpl,i)[1]:
	#return op(fplalt3(fpl,i))[1]:
fi:

wfpl := wmapfpl(fpl)[2*n-i+1]:
#nfpl := op(fplalt3(wfpl,2*n))[1]:
nfpl := fplalt3(wfpl,2*n)[1]:
return wmapfpl(nfpl)[i+1]:

end:


#fplalt2(fpl,i) lists the alt paths got from starting point i in the fpl
fplalt2 := proc(fpl,i) local j,k,B,W,W2,W3,S3,S1,S2,S,T,T2,n,ipt,ip1,conn,iconn,cnt,d:

n := nops(fpl):
conn := fplconn(fpl):

iconn := ei(conn,i):

if iconn=conn then
	return [[fpl,0]]:
#	return [fpl]:
fi:

S := {seq(seq({[j,k],[j+1,k]},j=0..n),k=1..n)}
union {seq(seq({[j,k],[j,k+1]},k=0..n),j=1..n)}:

B := fpltograph(fpl)[2]:
W := S minus B:


if i < 2*n then
	ip1:=i+1:
else
	ip1:=1:
fi:

S := distpaths2(findpath(fpl,i),findpath(fpl,ip1)):
S1 := []:
for j from 1 to nops(S) do
	if abs(S[j][1][1][1]-S[j][2][1][1]) + abs(S[j][1][1][2]-S[j][2][1][2]) = 1 then
		S1 := [op(S1),[S[j][1][2],S[j][1][1],S[j][2][1],S[j][2][2]]]:
	fi:
	if abs(S[j][1][1][1]-S[j][2][2][1]) + abs(S[j][1][1][2]-S[j][2][2][2]) = 1 then
		S1 := [op(S1),[S[j][1][2],S[j][1][1],S[j][2][2],S[j][2][1]]]:
	fi:
	if abs(S[j][1][2][1]-S[j][2][1][1]) + abs(S[j][1][2][2]-S[j][2][1][2]) = 1 then
		S1 := [op(S1),[S[j][1][1],S[j][1][2],S[j][2][1],S[j][2][2]]]:
	fi:
	if abs(S[j][1][2][1]-S[j][2][2][1]) + abs(S[j][1][2][2]-S[j][2][2][2]) = 1 then
		S1 := [op(S1),[S[j][1][1],S[j][1][2],S[j][2][2],S[j][2][1]]]:
	fi:
od:
#return S1:

T := []:
W2 := []:
cnt := 0:
for j from 1 to nops(S1) while cnt=0 do

	S2 := [appath(fpl,S1[j])]:
	S3 := aps1(B,W,S1[j],S1[j][1]):
	for k from 1 to nops(S2) do
		if fplconn(S2[k]) = iconn and looparea(S3[k]) <= 2*n-1 then
			T := [op(T),S2[k]]:
			W2 := [op(W2),S3[k]]:
			cnt := 1:
		fi:
	od:
	if nops(T) > 1 then
		d := min(seq(nops(W2[k]),k=1..nops(W2))):
		W3 := []:
		T2 := []:
		for k from 1 to nops(W2) do
			if nops(W2[k]) = d then
				W3 := [op(W3),W2[k]]:
				T2 := [op(T2),T[k]]:
			fi:
		od:
		T:=T2:
		W2:=W3:
		if nops(W3)>1 then
			print(fpl,i,W3):
		fi:
	fi:

od:

if nops(T)=0 then
	print(`None`,fpl,i,W3):
fi:

#T := [seq([T[k],nops(W2[k])],k=1..nops(T2))]:
T := [seq([T[k],looparea(W2[k])],k=1..nops(T))]:

return T:

end:

#fplalt3(fpl,i) lists the alt paths got from starting point i in the fpl
fplalt3 := proc(fpl,i) local j,k,B,W,W2,W3,S3,S1,S2,S,T,T2,n,ipt,ip1,conn,iconn,cnt,d:

n := nops(fpl):
conn := fplconn(fpl):
iconn := ei(conn,i):

if iconn=conn then
	return [fpl]:
fi:

S := {seq(seq({[j,k],[j+1,k]},j=0..n),k=1..n)}
union {seq(seq({[j,k],[j,k+1]},k=0..n),j=1..n)}:

B := fpltograph(fpl)[2]:
W := S minus B:


if i < 2*n then
	ip1:=i+1:
else
	ip1:=1:
fi:

S := distpaths2(findpath(fpl,i),findpath(fpl,ip1)):
S1 := []:
for j from 1 to nops(S) do
	if abs(S[j][1][1][1]-S[j][2][1][1]) + abs(S[j][1][1][2]-S[j][2][1][2]) = 1 then
		S1 := [op(S1),[S[j][1][2],S[j][1][1],S[j][2][1],S[j][2][2]]]:
	fi:
	if abs(S[j][1][1][1]-S[j][2][2][1]) + abs(S[j][1][1][2]-S[j][2][2][2]) = 1 then
		S1 := [op(S1),[S[j][1][2],S[j][1][1],S[j][2][2],S[j][2][1]]]:
	fi:
	if abs(S[j][1][2][1]-S[j][2][1][1]) + abs(S[j][1][2][2]-S[j][2][1][2]) = 1 then
		S1 := [op(S1),[S[j][1][1],S[j][1][2],S[j][2][1],S[j][2][2]]]:
	fi:
	if abs(S[j][1][2][1]-S[j][2][2][1]) + abs(S[j][1][2][2]-S[j][2][2][2]) = 1 then
		S1 := [op(S1),[S[j][1][1],S[j][1][2],S[j][2][2],S[j][2][1]]]:
	fi:
od:

S2 := [seq([appath(fpl,S1[j])],j=1..nops(S1))]:
S3 := [seq(aps1(B,W,S1[j],S1[j][1]),j=1..nops(S1))]:

#return S2,S3:

T := []:
W2 := []:
cnt := 0:
for j from 1 to nops(S2) do
	for k from 1 to nops(S2[j]) do
		if fplconn(S2[j][k]) = iconn then
			T := [op(T),S2[j][k]]:
			W2 := [op(W2),S3[j][k]]:
		fi:
	od:
od:

#return T,W2:


if nops(T) > 1 then
	d := min(seq(nops(W2[k]),k=1..nops(W2))):
	W3 := []:
	T2 := []:
	for k from 1 to nops(W2) do
		if nops(W2[k]) = d then
			W3 := [op(W3),W2[k]]:
			T2 := [op(T2),T[k]]:
		fi:
	od:
	T := T2:
fi:

return T:
#return [T[1]]:

end:

#fplalt33(fpl,i) lists the alt paths got from starting point i in the fpl
fplalt33 := proc(fpl,i) local j,k,B,W,W2,W3,S3,S1,S2,S,T,T2,n,ipt,ip1,conn,iconn,cnt,d:

n := nops(fpl):
conn := fplconn(fpl):
iconn := ei(conn,i):

if iconn=conn then
	return [fpl,i]:
fi:

S := {seq(seq({[j,k],[j+1,k]},j=0..n),k=1..n)}
union {seq(seq({[j,k],[j,k+1]},k=0..n),j=1..n)}:

B := fpltograph(fpl)[2]:
W := S minus B:


if i < 2*n then
	ip1:=i+1:
else
	ip1:=1:
fi:

S := distpaths2(findpath(fpl,i),findpath(fpl,ip1)):
S1 := []:
for j from 1 to nops(S) do
	if abs(S[j][1][1][1]-S[j][2][1][1]) + abs(S[j][1][1][2]-S[j][2][1][2]) = 1 then
		S1 := [op(S1),[S[j][1][2],S[j][1][1],S[j][2][1],S[j][2][2]]]:
	fi:
	if abs(S[j][1][1][1]-S[j][2][2][1]) + abs(S[j][1][1][2]-S[j][2][2][2]) = 1 then
		S1 := [op(S1),[S[j][1][2],S[j][1][1],S[j][2][2],S[j][2][1]]]:
	fi:
	if abs(S[j][1][2][1]-S[j][2][1][1]) + abs(S[j][1][2][2]-S[j][2][1][2]) = 1 then
		S1 := [op(S1),[S[j][1][1],S[j][1][2],S[j][2][1],S[j][2][2]]]:
	fi:
	if abs(S[j][1][2][1]-S[j][2][2][1]) + abs(S[j][1][2][2]-S[j][2][2][2]) = 1 then
		S1 := [op(S1),[S[j][1][1],S[j][1][2],S[j][2][2],S[j][2][1]]]:
	fi:
od:

S2 := [seq([appath(fpl,S1[j])],j=1..nops(S1))]:
S3 := [seq(aps1(B,W,S1[j],S1[j][1]),j=1..nops(S1))]:

#return S2,S3:

T := []:
W2 := []:
cnt := 0:
for j from 1 to nops(S2) do
	for k from 1 to nops(S2[j]) do
		if fplconn(S2[j][k]) = iconn then
			T := [op(T),S2[j][k]]:
			W2 := [op(W2),S3[j][k]]:
		fi:
	od:
od:

#return T,W2:


if nops(T) > 1 then
	d := min(seq(nops(W2[k]),k=1..nops(W2))):
	W3 := []:
	T2 := []:
	for k from 1 to nops(W2) do
		if nops(W2[k]) = d then
			W3 := [op(W3),W2[k]]:
			T2 := [op(T2),T[k]]:
		fi:
	od:
	T := T2:
else
	W3 := W2:
fi:

if nops(W3[1]) mod 4 = 1 then
	if (i+(nops(W3[1])-1)/4) > 2*n then
		cnt := (i+(nops(W3[1])-1)/4) - 2*n:
	else
		cnt := (i+(nops(W3[1])-1)/4):
	fi:
elif nops(W3[1]) mod 4 = 3 then
	if (i+(nops(W3[1])-3)/4) > 2*n then
		cnt := (i+(nops(W3[1])-3)/4) -2*n: 
	else
		cnt := (i+(nops(W3[1])-3)/4):
	fi:
else 
	cnt := 0:
fi:

return [T[1],cnt]:

end:


#bumpfpl(fpl) chooses the new fpl according the bumping procedure with minimal displacement
bumpfpl := proc(fpl) local B1,B2,n,conn,iconn,i,j,k,f2,min1,S,T:

n := nops(fpl):
conn := fplconn(fpl):

T := []:

for i from 1 to 2*n do

	iconn := ei(conn,i):
	if iconn=conn then
		T := [op(T),fpl]:
	else
		f2 := npitofpl(iconn):
		B1 := fpltograph(fpl)[2]:
		min1 := nops(B1):
		S := {}:

		for j from 1 to nops(f2) do
			B2 := fpltograph(f2[j])[2]:
			if nops(B1 minus B2) < min1 and not evalb(f2[j] in {op(T)}) then
				S := {f2[j]}:
				min1 := nops(B1 minus B2):
			elif nops(B1 minus B2) = min1 and not evalb(f2[j] in {op(T)}) then
				S := S union {f2[j]}:
			fi:
		od:

		if nops(S) > 1 then 
			print(fpl,i,`More`,nops(S)):
		elif nops(S) = 0 then 
			print(fpl,i,`None`):
		fi:
		T := [op(T),op(S)]:
	fi:
od:

return T:

end:

#wmapalt(fpl) checks if fplalt(fpl,i) rotated returns the same fpl as fplalt(rotated fpl, i+1)
wmapalt := proc(fpl) local i,j,wfpl,f2,w2,n:

n := nops(fpl):
wfpl := wmapfpl(fpl)[2]:
f2 := [seq(wmapfpl(op(fplalt3(fpl,i)))[2],i=1..2*n)]:
w2 := [seq(op(fplalt3(wfpl,i)),i=2..2*n),op(fplalt3(wfpl,1))]:
j := 0:
for i from 1 to 2*n do
	if f2[i] <> w2[i] then
		print(`Not equal`,fpl,i,wfpl,i+1,f2[i],w2[i]):
		j := 1:
	fi:
od:

if j=0 then
	print(`All fine`):
fi:

end:

#brfr(n) tries the brute force algorithm for determining the correct action of e_{2n}
brfr1 := proc(n) local i,j,fs,nfs,asms,conns,iconn,S,S2,Sbar,A,T:

conns := numchords(1,2*n):
S := []:
S2 := []:

for i from 1 to nops(conns) do
	iconn := ei(conns[i],2*n):
	fs := npitofpl(conns[i]):
	nfs := npitofpl(iconn):
	if iconn=conns[i] then
		S := [op(S),seq([fs[j],fs[j]],j=1..nops(fs))]:
	elif nops(nfs)=1 then
		S := [op(S),seq([fs[j],nfs[1]],j=1..nops(fs))]:
	else
		S2 := [op(S2),choose1(fs,nfs)]:
	fi:
od:

#return S,S2:
S2 := choose2(S2):

A := {}:

for i from 1 to nops(S2) do
	Sbar := [op(S),op(S2[i])]:
	T := convert([seq(op(wmapfpl(Sbar[j][2])),j=1..nops(Sbar))],multiset):
	if {seq(T[j][2],j=1..nops(T))} = {2*n} then
		A := A union {Sbar}:
	fi:
od:

A:

end:

#choose1(L1,L2) returns a list of all possible ways of matching them
choose1 := proc(L1,L2) local i,j,m,n,Ve2,Ve,S:

S := []:
n:=nops(L1):
m:=nops(L2):

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

for i from 1 to nops(Ve2) do
	if nops(Ve2[i]) < n then
		Ve := Ve union {[op(Ve2[i]),0$(n-nops(Ve2[i]))]}:
	else
		Ve := Ve union {Ve2[i]}:
	fi:
od:

for i from 1 to nops(Ve) do
	S := [op(S),[seq([L1[j],L2[Ve[i][j]+1]],j=1..nops(L1))]]:
od:

return S:

end:

choose2 := proc(S) local i,j,L,m,n,Ve2,Ve,cnt:

L := [seq(nops(S[i]),i=1..nops(S))]:
n:= nops(L):
m := max(op(L)):

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

for i from 1 to nops(Ve2) do
	if nops(Ve2[i]) < n then
		Ve := Ve union {[op(Ve2[i]),0$(n-nops(Ve2[i]))]}:
	else
		Ve := Ve union {Ve2[i]}:
	fi:
od:

Ve2 := Ve:
Ve:= {}:
for i from 1 to nops(Ve2) do
	cnt := 0:
	for j from 1 to nops(L) do
		if Ve2[i][j]+1 > L[j] then
			cnt := 1:
		fi:
	od:
	if cnt=0 then
		Ve := Ve union {Ve2[i]}:
	fi:
od:

return [seq([seq(op(S[j][Ve[i][j]+1]),j=1..nops(Ve[i]))],i=1..nops(Ve))]:

end:


######################################Procedures about counting FPLs by looking at Lattice Paths#########


#roads(stpart,ept,bound) given the starting part of the path, the ending point and the boundary returns the set of all possible non-intersecting paths within the boundary
roads := proc(stpart,ept,bound) local i,j,p1,p2,rds,spt:
option remember:

if stpart[nops(stpart)] = ept then
	return {stpart}:
fi:

if evalb(stpart[nops(stpart)] in bound) then
	return {}:
fi:

spt := stpart[nops(stpart)]:

p1 := {[spt[1]+1,spt[2]],[spt[1]-1,spt[2]],[spt[1],spt[2]+1],[spt[1],spt[2]-1]}:
p1 := p1 minus convert(stpart,set):

p2 := {}:

for i from 1 to nops(p1) do
	p2 := p2 union roads([op(stpart),p1[i]],ept,bound):
od:

return p2:

end:


#pitoppl(conn) given a connectivity finds all possible partially packed loops
pitoppl := proc(conn) local i,j,bdry,dim,rdconn,paths,lnodes:

dim := nops(conn):

bdry := {seq([0,i],i=0..dim+1),seq([i,0],i=0..dim+1),seq([dim+1,i],i=0..dim+1),seq([i,dim+1],i=0..dim+1)}:
bdry := bdry union {seq([1,2*i],i=1..floor(dim/2))} union {seq([dim,dim+1-2*i],i=1..floor(dim/2))} union {seq([2*ceil(dim/2)+1-2*i,1],i=1..ceil(dim/2))} union {seq([2*i+dim-2*ceil(dim/2),dim],i=1..ceil(dim/2))}:

rdconn := subs({seq(i=[1,2*i],i=1..floor(dim/2)),
		seq(floor(dim/2)+i=[2*i+dim-2*ceil(dim/2),dim],i=1..ceil(dim/2)),
		seq(dim+i=[dim,dim+1-2*i],i=1..floor(dim/2)),
		seq(dim+floor(dim/2)+i=[2*ceil(dim/2)+1-2*i,1],i=1..ceil(dim/2))},conn):

paths := parentroad(rdconn,bdry):


return paths:

end:

#pitofpl(conn) given a connectivity finds all possible fpls
pitofpl := proc(conn) local i,j,bdry,dim,rdconn,paths,paths2,lnodes,allnodes,lpnodes:

dim := nops(conn):

#cc :=subs({seq([0,2*i]=i,i=1..floor(dim/2))},cc):
#cc :=subs({seq([2*i+dim-2*ceil(dim/2),dim+1]=floor(dim/2)+i,i=1..ceil(dim/2))},cc):
#cc :=subs({seq([dim+1,dim+1-2*i]=dim+i,i=1..floor(dim/2))},cc):
#cc :=subs({seq([2*ceil(dim/2)+1-2*i,0]=dim+floor(dim/2)+i,i=1..ceil(dim/2))},cc):

bdry := {seq([0,i],i=0..dim+1),seq([i,0],i=0..dim+1),seq([dim+1,i],i=0..dim+1),seq([i,dim+1],i=0..dim+1)}:
bdry := bdry union {seq([1,2*i],i=1..floor(dim/2))} union {seq([dim,dim+1-2*i],i=1..floor(dim/2))} union {seq([2*ceil(dim/2)+1-2*i,1],i=1..ceil(dim/2))} union {seq([2*i+dim-2*ceil(dim/2),dim],i=1..ceil(dim/2))}:

rdconn := subs({seq(i=[1,2*i],i=1..floor(dim/2)),
		seq(floor(dim/2)+i=[2*i+dim-2*ceil(dim/2),dim],i=1..ceil(dim/2)),
		seq(dim+i=[dim,dim+1-2*i],i=1..floor(dim/2)),
		seq(dim+floor(dim/2)+i=[2*ceil(dim/2)+1-2*i,1],i=1..ceil(dim/2))},conn):

paths := parentroad(rdconn,bdry):
paths2 := {}:
allnodes := {seq(seq([i,j],i=1..dim),j=1..dim)}:

for i from 1 to nops(paths) do
	lnodes := allnodes minus {seq(seq(paths[i][j][k],k=1..nops(paths[i][j])),j=1..nops(paths[i]))}:
	if nops(lnodes) = 0 then
		paths2 := paths2 union {paths[i]}:
	else
		lpnodes := fillloop(lnodes):
		if nops(lpnodes)>0 and lpnodes <> {{}} then
			paths2 := paths2 union {seq({op(paths[i]),lpnodes[j]},j=1..nops(lpnodes))}:
		fi:
	fi:
od:

return paths2:

end:


childroad := proc(spt,ept,bdry) local bdr:
option remember:

bdr := bdry minus {spt}:

return roads([spt],ept,bdr):

end:

parentroad := proc(rdconn,bdry) local i,j,croad,proad,p2road:
option remember:

if rdconn={} then
	return {{}}:
fi:

proad := {}:

croad := childroad(rdconn[1][1],rdconn[1][2],bdry):

for i from 1 to nops(croad) do
	p2road := parentroad(rdconn minus {rdconn[1]},bdry union {op(croad[i])}):
	proad := proad union {seq({op(p2road[j]),croad[i]},j=1..nops(p2road))}:
od:

return proad:

end:

#fillloops(sites) finds all possible sets of loops which fill all sites
fillloops := proc(sites) local i,j,k,s2,cnt,lps:

s2 := {}:
for i from 1 to nops(sites) do
	cnt:=0:
	for j from 1 to nops(s2) do
		for k from 1 to nops(s2[j]) do
			if sites[i]=[s2[j][k][1]+1,s2[j][k][2]] or sites[i]=[s2[j][k][1]-1,s2[j][k][2]] or sites[i]=[s2[j][k][1],s2[j][k][2]+1] or sites[i]=[s2[j][k][1],s2[j][k][2]-1] then
				s2 := {op(1..j-1,s2),{op(1..nops(s2[j]),s2[j]),sites[i]},op(j+1..nops(s2),s2)}:
				cnt:=1:
			fi:
		od:
	od:
	if cnt=0 then
		s2 := {op(s2),{sites[i]}}:
	fi:
od:

lps := {seq(fillloop(s2[i]),i=1..nops(s2))}:

return s2:

end:

#fillloop(sites) finds all loops which fills all sites and come back home
fillloop := proc(sites) local lps,lps2,lps3,i,j:
option remember:

if nops(sites) < 4 then
	return {{}}:
fi:

lps2 := {}:

for i from 4 to nops(sites) do
	lps := distn(sites,sites[1],i):
	for j from 1 to nops(lps) do
		if lps[j][nops(lps[j])] = sites[1] then
			lps2 := lps2 union {lps[j]}:
		fi:
	od:
od:

lps := lps2:
lps2 := {}:

for i from 1 to nops(lps) do
        if not evalb([lps[i][1],seq(lps[i][nops(lps[i])+1-j],j=2..nops(lps[i]))] in lps2) then
                lps2 := lps2 union {lps[i]}:
        fi:
od:

#return lps2:

lps := {}:

for i from 1 to nops(lps2) do
	lps3 := fillloop(sites minus {op(lps2[i])}):
	lps := lps union {seq({lps2[i],op(lps3[j])},j=1..nops(lps3))}:
od:

lps2 := {}:
for i from 1 to nops(lps) do
	if {seq(seq(lps[i][j][k],k=1..nops(lps[i][j])),j=1..nops(lps[i]))} = sites then
		lps2 := lps2 union {lps[i]}:
	fi:
od:

return lps2:

end:

fllp := proc(sites,pathsofar) local i,j,n,lastpt,S:
option remember:

n := nops(pathsofar):
lastpt := pathsofar[n]:

S := {}:

if nops(sites)=0 then
	if ([lastpt[1]+1,lastpt[2]]=pathsofar[1] or [lastpt[1]-1,lastpt[2]]=pathsofar[1] or [lastpt[1],lastpt[2]+1]=pathsofar[1] or [lastpt[1],lastpt[2]-1]=pathsofar[1]) and nops(pathsofar)>3 then
		return {pathsofar}:
	else
		return {}:
	fi:
fi:

for i from 1 to nops(sites) do
	if [lastpt[1]+1,lastpt[2]]=sites[i] or [lastpt[1]-1,lastpt[2]]=sites[i] or [lastpt[1],lastpt[2]+1]=sites[i] or [lastpt[1],lastpt[2]-1]=sites[i] then
		S := S union {[op(pathsofar),sites[i]]}:
	fi:
od:

return {seq(op(fllp(sites minus {S[i][nops(S[i])]},S[i])),i=1..nops(S))}:

end:


dist1:=proc(sites,v) local i,S:

S:={}:

for i from 1 to nops(sites) do
	if [v[1]+1,v[2]]=sites[i] or [v[1]-1,v[2]]=sites[i] or 
	[v[1],v[2]+1]=sites[i] or [v[1],v[2]-1]=sites[i] then 
		S := S union {sites[i]}:
	fi:
od:

return S:

end:

distn:=proc(sites,v,n) local i,j,k,S,T,S2:

S := {[v]}:

for i from 1 to n do
	S2 := {}:
	for j from 1 to nops(S) do
		T := dist1(sites,S[j][nops(S[j])]):
		S2 := S2 union {seq([op(S[j]),T[k]],k=1..nops(T))}:
	od:
	S := S2:
od:

S2 := {}:

for i from 1 to nops(S) do
	if nops({op(1..nops(S[i])-1,S[i])}) = nops([op(1..nops(S[i])-1,S[i])]) and nops({op(2..nops(S[i]),S[i])}) = nops([op(2..nops(S[i]),S[i])]) then
		S2 := S2 union {S[i]}:
	fi:
od:

return S2:

end:

##################Procedures about chordal relationships in link patterns#######
#dset(pat) returns the distance set of a given pattern
dset := proc(pat) local i,ans:

ans := []:

for i from 1 to nops(pat) do
	if abs(pat[i][1]-pat[i][2]) < nops(pat) then
		ans:= [op(ans),abs(pat[i][1]-pat[i][2])]:
	else
		ans:= [op(ans),2*nops(pat)-abs(pat[i][1]-pat[i][2])]:
	fi:
od:

return convert(ans,multiset):

end:



#part2(S) given a set S, returns the set of all possible 2-pairings of the set
part2 := proc(S) local i,S2,ans:
option remember:

if nops(S) mod 2 = 1 then
	ERROR(`The set must have an even number of elements.`):
fi:

if nops(S) = 2 then
	return {{{S[1],S[2]}}}:
fi:

ans := {}:

S2 := combinat[choose](S,2):

for i from 1 to nops(S2) do
	ans := ans union {seq({S2[i],op(part2(S minus S2[i])[j])},j=1..nops(part2(S minus S2[i])))}:
od:

return ans:

end:

#labchords(n) depicts the non-intersecting chords on the labelled set {1,...,2n}
labchords := proc(n) local S,i,j,k,C,C2,cij,cik:

S := {seq(i,i=1..2*n)}:

C := part2(S):

C2 := C:

for i from 1 to nops(C) do
	for j from 1 to nops(C[i]) do
		for k from 1 to nops(C[i]) do

			cij := sort([op(C[i][j])]):
			cik := sort([op(C[i][k])]):

			if cik[1] < cij[1] then
				if cik[2] < cij[2] and cik[2] > cij[1] then
					C2 := C2 minus {C[i]}:
				fi:
			elif cik[1] < cij[2] then
				if cik[2] > cij[2] then
					C2 := C2 minus {C[i]}:
				fi:
			fi:
		od:
	od:
od:

return C2:

end:

#add1(S) adds one to each element in the given bipartition
add1 := proc(S) local n,i,j,S2:

n := nops(S):

S2 := [[0$2]$n]:

for i from 1 to n do
	for j from 1 to 2 do
		if S[i][j]=2*n then
			S2[i][j] := 1:
		else
			S2[i][j] := S[i][j]+1:
		fi:
	od:
od:

return {seq({op(S2[i])},i=1..n)}:

end:			



#unlabchords(n) returns the distance-set of chords on unlabelled points along with the number of times it appears in labchords(n)
unlabchords := proc(n) local S,C,C2,D,i,j:

C := labchords(n):

C2 := [seq(dset(C[i]),i=1..nops(C))]:

S := {op(C2)}:

D := [0$nops(S)]:

for i from 1 to nops(C2) do
	for j from 1 to nops(S) do
		if C2[i] = S[j] then
			D[j] := D[j]+1:
		fi:
	od:
od:

return seq({S[i],D[i]},i=1..nops(S)):


end:


#numchords(n1,n2) returns the non-intersecting chords of the points {n1,...,n2} on the circle
numchords := proc(n1,n2) local k,i,j,S,s1,s2:
option remember:

if n2 <= n1 then
	return {}:
fi:

if n2-n1=1 then
	return {{{n1,n2}}}:
fi:

#S := {seq(seq(seq({{n1,2*k},numchords(n1+1,2*k-1)[i],numchords(2*k+1,n2)[j]},i=1..nops(numchords(n1+1,2*k-1))),j=1..nops(numchords(2*k+1,n2))),k=floor(n1/2)+1..floor(n2/2))}:

S := {}:

for k from n1+1 to n2 by 2 do

	s1:=numchords(n1+1,k-1):
	s2:=numchords(k+1,n2):

	if nops(s1) = 0 then
		if nops(s2)=0 then
			S := S union {{n1,k}}:
		else
			S := S union {seq({{n1,k},op(s2[i])},i=1..nops(s2))}:
		fi:
	else
		if nops(s2)=0 then
			S := S union {seq({{n1,k},op(s1[i])},i=1..nops(s1))}:
		else
			S := S union {seq(seq({{n1,k},op(s1[i]),op(s2[j])},i=1..nops(s1)),j=1..nops(s2))}:
		fi:
	fi:
od:

return S:

end:

#ineqchords(n) returns the inequivalent chords on 2n points in a circle
ineqchords := proc(n) local i,j,k,S,S2,cnt:

S := numchords(1,2*n):

S2 := {}:

for i from 1 to nops(S) do
	cnt:=0:
	for j from 1 to 2*n-1 do
		if evalb(subs({seq(k=k+j,k=1..2*n-j),seq(k=k+j-2*n,k=2*n-j+1..2*n)},S[i]) in S2) then
			cnt:=1:
		fi:
	od:
	if cnt=0 then
		S2 := S2 union {S[i]}:
	fi:
od:

return S2:


end:



######################################Procedures about Height Functions#########
#asmtohtfn(asm) given an asm returns the corresponding height function
asmtohtfn := proc(asm) local i,j,k,l,htfn,n:

n := nops(convert(asm,listlist)):

htfn := [[0$(n+1)]$(n+1)]:

for i from 1 to n+1 do
	for j from 1 to n+1 do
		htfn[i][j] := i+j-2-2*add(add(asm[k-1,l-1],l=2..j),k=2..i):
	od:
od:

return htfn:

end:

#htfntoasm(htfn) given a height function returns the corresponding asm
htfntoasm := proc(htfn) local i,j,k,l,asm,n,cor:

n:= nops(htfn):

cor := [seq([seq((i+j-2-htfn[i][j])/2,j=1..n)],i=1..n)]:

asm := [seq([seq(cor[i][j]+cor[i-1][j-1]-cor[i][j-1]-cor[i-1][j],j=2..n)],i=2..n)]:

return asm:

end:

#wmap(htfn) acts on the height function htfn by the map G defined in Wieland's paper and returns the new height function
wmap := proc(htfn) local i,j,h2,n:

n := nops(htfn):
h2 := htfn:

for i from 2 to n-1 do
	for j from 2 to n-1 do
		if (i+j) mod 2 =0 then
			if h2[i][j-1]=h2[i][j+1] and h2[i][j-1]=h2[i+1][j] and h2[i][j-1]=h2[i-1][j] then
				if h2[i][j] < h2[i][j-1] then
					h2[i][j] := h2[i][j]+2:
				else
					h2[i][j] := h2[i][j]-2:
				fi:
			fi:
		fi:
	od:
od:

for i from 2 to n-1 do
	for j from 2 to n-1 do
		if (i+j) mod 2 = 1 then
			if h2[i][j-1]=h2[i][j+1] and h2[i][j-1]=h2[i+1][j] and h2[i][j-1]=h2[i-1][j] then
				if h2[i][j] < h2[i][j-1] then
					h2[i][j] := h2[i][j]+2:
				else
					h2[i][j] := h2[i][j]-2:
				fi:
			fi:
		fi:
	od:
od:


return h2:

end:

#wmapfpl(fpl) given an fpl finds all other fpls which are related to it by rotation
wmapfpl := proc(fpl) local htfn,i,j,S:

htfn := wmap(asmtohtfn(fpltoasm(fpl))):
S := [fpl]:
#printfpl(fpl):

#for i from 1 while asmtofpl(htfntoasm(htfn)) <> fpl do
for i from 1 to 2*nops(fpl)-1 do
	#print(i):
	#printfpl(asmtofpl(htfntoasm(htfn))):
	S := [op(S),asmtofpl(htfntoasm(htfn))]:
	htfn := wmap(htfn):
od:

#printfpl(asmtofpl(htfntoasm(wmap(htfn)))):

return S:

end:




######################################Procedures from CLIQUE written for Math640 in 2006 by Doron Zeilberger#########

#Neig1(G,v): the set  of neighbors of vertex v in G
Neig1:=proc(G,v) local S,e,V,E:
S:={}: V:=G[1]: E:=G[2]:

for e in E do
 if member(v,e) then
  S:=S union  {(e minus {v})[1]} :
 fi:
od:

S:
end:


#Neig(G,S): the set  of all neighbors of set of vertices S
Neig:=proc(G,S) local V,s: 

{seq(op(Neig1(G,s)), s in S)}:
end:



##CC(G,v): Inputs a simple graph G and a vertex v
#outputs the set of vertices that can be reached from v
CC:=proc(G,v) local V,E,S1,S2:

V:=G[1]: E:=G[2]:
S1:={v}:

S2:= {v} union Neig(G,S1):

while S1<>S2 do
 S1:=S2:
 S2:=S2 union Neig(G,S2):
od:
S2:
end:


#CCD(G): inputs a graph and outputs its set of connected
#components
CCD:=proc(G) local V,E,S,cc,v:
V:=G[1]: E:=G[2]:

S:={}:

while V<>{} do
 v:=V[1]:

 cc:=CC(G,v):
 S:= S union {cc}:
 V:= V minus cc:
od:

S:
end:



##################Procedures from ROBBINS by Doron Zeilberger#########
ASM:=proc(k)
local i,gu,asm:
 
gu:=GOGset(k,k):
 
print(`There are`, nops(gu),`Alternating Sign Matrices of size`,k):
print(`Here they all are:`):
 
for i from 1 to nops(gu) do
  asm:=GOGTOASM(op(i,gu)):
 print(op(asm)):
od:
 
 
end:

 
#GOGa(k,n,a)  gives  the   set of  k by n Gog-trapezoids  such  that
#the rightmost border is the  vector a
 
GOGa:=proc(k,n,a)
local pip,kvu,firow,b,mu,gu,i,j,l,trap,trap1:
 
if  not k>=1 or not n>=k  or  not nops(a)=k  then
  ERROR(`Improper intput`):
fi:
 
 
if n=k and k=1 then
   if not op(1,a)=1 then
     RETURN({}):
   else
     RETURN({[[1]]}):
   fi:
fi:
 
if n=k then
   if not op(1,a)=k then
      ERROR(`Wrong input`):
   fi:
   mu:=GOGa(k-1,k,[op(2..k,a)]):
    gu:={}:
 
    for i from 1 to nops(mu) do
     trap1:=op(i,mu):
     firow:=op(1,trap1):
    
     gu:=gu union {[[op(firow),k],op(2..k,trap1)]}:
        od:
 RETURN(gu):
fi:
 
gu:={}:    
 
kvu:=Tkn(k,n,a):
 
 
 for pip from 1 to nops(kvu) do
   b:=op(pip,kvu):
 
  mu:=GOGa(k,n-1,b):
 
   for j from 1 to nops(mu) do
    trap:=op(j,mu):
    trap1:=[op(1..n-k,trap)]:
 
     for l from n-k+1 to n-1 do
       trap1:=[op(trap1),[op(op(l,trap)),op(l-(n-k),a)]]:
     od:
      
      trap1:=[op(trap1),[op(k,a)]]:
     gu:=gu union {trap1}:
 
   od:
 
 
od:
 
gu:
 
end:
 
GOG:=proc(k,n)
local i,gu,lu:
 
lu:=LOGOG(k,n):
gu:={}:
 
for i from 1 to nops(lu) do
 gu:=gu union GOGa(k,n,op(i,lu)):
od:
 
gu:
 
print(`The number of Gog Trapezoids with k=`,k,`and n=`,n,`equals`,nops(gu)):
print(`Here they all are`):
 
for i from 1 to nops(gu) do
yafe(op(i,gu)):
lprint(``):
od:
 
gu:
 
end:
 
 
GOGset:=proc(k,n)
local i,gu,lu:
 
lu:=LOGOG(k,n):
gu:={}:
 
for i from 1 to nops(lu) do
 gu:=gu union GOGa(k,n,op(i,lu)):
od:
 
gu:
 
gu:
 
end:
 
GOGTOASM:=proc(mt)
local k,mat,mat1,ro,i,j:
 
k:=nops(mt):
 
mat:=array(1..k,1..k):
mat1:=array(1..k,1..k):
 
for i from 1 to k do
 for j from 1 to k do
   mat[i,j]:=0:
 od:
od:
 
 
  for i from 1 to k do
   ro:=op(i,mt):
    for j from 1 to nops(ro) do
      mat[i,op(j,ro)]:=1:
    od:
  od:
 
for i from 1 to k-1 do
 for j from 1 to k do
  mat1[i,j]:=mat[i,j]-mat[i+1,j]:
 od:
 
 for j from 1 to k do
  mat1[k,j]:=mat[k,j]:
 od:
 
od:
 
mat1:
 
end:

 
DECSEQ:=proc(k,n)
local gu,gu1,i1:
option remember:
 
if n=1 and  k=1 then
 RETURN({[1]}):
fi:
 
if k=0  and  n>=1 then
    RETURN({[]}):
fi:
 
if n<1 or  k<1  then
   RETURN({}):
fi:
 
gu:=DECSEQ(k,n-1):
gu1:=DECSEQ(k-1,n):
 
for i1 from 1   to  nops(gu1)  do
  gu:=gu  union   {[n,op(op(i1,gu1))]}:
od:
 
gu:
 
end:
 
 
DECSEQ0:=proc(k,n)
local gu,gu1,i1:
option remember:
 
if n=0 and  k=1 then
 RETURN({[0]}):
fi:
 
if k=0  and  n>=1 then
    RETURN({[]}):
fi:
 
if n<0 or  k<1  then
   RETURN({}):
fi:
 
gu:=DECSEQ0(k,n-1):
gu1:=DECSEQ0(k-1,n):
 
for i1 from 1   to  nops(gu1)  do
  gu:=gu  union   {[n,op(op(i1,gu1))]}:
od:
 
gu:
 
end:
 
 
LOGOG:=proc(k,n)
local gu,gu1,i,i1,vec,nakh:
 
if not k>=1  or  not n>=k then
  RETURN({}):
fi:
 
gu1:=DECSEQ(k,n):
 
gu:={}:
 
 
for i from 1 to nops(gu1) do
 vec:=op(i,gu1):
  nakh:=1:
    for  i1  from 1 to k  do
       if  not op(i1,vec)>=k-i1+1 then
           nakh:=0:
           exit:
       fi:
    od:
 
      if nakh=1 then
          gu:=gu union {vec}:
        fi:
od:
 
 
gu:
 
end:
 
 
ELOGOG:=proc(k,n)
local gu,gu1,i,i1,vec,nakh:
 
if not k>=1  or  not n>=k then
  RETURN({}):
fi:
 
gu1:=DECSEQ0(k,n+1):
 
gu:={}:
 
 
for i from 1 to nops(gu1) do
 vec:=op(i,gu1):
  nakh:=1:
 
    if op(1,vec)=n+1 and op(2,vec)>n then
      nakh:=0:
    fi:
 
    if op(1,vec)=n and n=k and op(2,vec)=k then
      nakh:=0:
    fi:
 
    for  i1  from 1 to k  do
       if  not op(i1,vec)>=k-i1 then
           nakh:=0:
           exit:
       fi:
    od:
 
      if nakh=1 then
          gu:=gu union {vec}:
        fi:
od:
 
 
gu:
 
end:
 
 
 
#Tkn gives the set T_k(n;a), where a=[a_1, ..., a_k] defined in the
#proof of 1.2.1.1
 
Tkn:=proc(k,n,a)
local nakh,i1,b,i,gu,mu:
 
mu:=DECSEQ(k,n-1):
gu:={}:
 
for i from 1 to nops(mu) do
  b:=op(i,mu):
 
nakh:=1:
 
 
 
  for i1 from 1 to 1 do
    if not ( k-i1+1<=op(i1,b) and op(i1,b)<=op(i1,a) )
        then
            nakh:=0:
    fi:
 
  od:
 
  for i1 from 2 to k do
    if not ( k-i1+1<=op(i1,b) and op(i1,b)<=min( op(i1,a),op(i1-1,a)-1 )   )
        then
            nakh:=0:
    fi:
 
  od:
 
 
if nakh=1 then
  gu:=gu union {b}:
fi:
 
od:
 
gu:
 
end:




############################################Rejected Algorithms##########################################







#exfplalt(fpl,i) lists the alt paths got from starting point i in the fpl
exfplalt := proc(fpl,i) local j,k,B,W,S1,S2,S,n,ipt,ip1,conn,iconn,ap1,ap2,stpart,sq,fs1,fs2:

n := nops(fpl):
conn := fplconn(fpl):
iconn := ei(conn,i):

if iconn=conn then
	return {fpl}:
fi:

S := {seq(seq({[j,k],[j+1,k]},j=0..n),k=1..n)}
union {seq(seq({[j,k],[j,k+1]},k=0..n),j=1..n)}:

B := fpltograph(fpl)[2]:

W := S minus B:

if i <= floor(n/2) then
	ipt := [1,2*i]:
elif i <= n then
	ipt := [2*i-n,n]:
elif i <= n + floor(n/2) then
	ipt := [n,3*n+1-2*i]:
else
	ipt := [4*n+1-2*i,1]:
fi:

if i < floor(n/2) then
	ip1 := [1,2*(i+1)]:
elif i < n then
	ip1 := [2*(i+1)-n,n]:
elif i < n + floor(n/2) then
	ip1 := [n,3*n+1-2*(i+1)]:
elif i < 2*n then
	ip1 := [4*n+1-2*(i+1),1]:
else
	ip1 := [1,2]:
fi:

ip1 := ipt:

if i < floor(n/2) then
	stpart := [ipt,[1,ipt[2]+1],[1,ipt[2]+2]]:
	sq := {ipt,[1,ipt[2]+1]}:
elif i = floor(n/2) and n mod 2 = 1 then
	stpart := [ipt,[1,ipt[2]+1],[2,ipt[2]+1]]:
	sq := {ipt}:
elif i = floor(n/2) and n mod 2 = 0 then
	stpart := [[1,ipt[2]-1],ipt,[2,ipt[2]]]:
	sq := {[1,ipt[2]-1]}:
	ip1 := [1,ipt[2]-1]:
elif i < n then
	stpart := [ipt,[ipt[1]+1,n],[ipt[1]+2,n]]:
	sq := {[ipt[1],n-1],[ipt[1]+1,n-1]}:
elif i = n  then
	stpart := [[n-1,n],ipt,[n,n-1]]:
	sq := {[n-1,n-1]}:
	ip1 := [n-1,n]:
elif i < n+floor(n/2) then
	stpart := [ipt,[n,ipt[2]-1],[n,ipt[2]-2]]:
	sq := {[n-1,ipt[2]-1],[n-1,ipt[2]-2]}:
elif i = n+floor(n/2) and n mod 2 = 0 then
	stpart := [[n,2],ipt,[n-1,1]]:
	sq := {[n-1,1]}:
	ip1 := [n,2]:
elif i = n+floor(n/2) and n mod 2 = 1 then
	stpart := [ipt,[n,1],[n-1,1]]:
	sq := {[n-1,1]}:
elif i < 2*n then
	stpart := [ipt,[ipt[1]-1,1],[ipt[1]-2,1]]:
	sq := {[ipt[1]-1,1],[ipt[1]-2,1]}:
elif i = 2*n  then
	stpart := [[2,1],ipt,[1,2]]:
	sq := {ipt}:
	ip1 := [2,1]:
fi:

#altpaths := aps1(B,W,stpart,ip1):
#if evalb(iconn in {seq(fplconn(newaltfpl(fpl,altpaths[j])),j=1..nops(altpaths))})  then
#	print(`Yes`):
#else
#	print(`No`):
#fi:
#return altpaths:

S1 := {}:
S2 := {}:

if nops(sq) = 1 then
	ap1 := apsq(fpl,sq[1]):
	fs1 := [seq(newaltfpl(fpl,ap1[j]),j=1..nops(ap1))]:
	for j from 1 to nops(fs1) do
		if fplconn(fs1[j])=iconn then
			S1 := S1 union {fs1[j]}:
		fi:
	od:
	if nops(S1) = 0 then
		print(`nothing`,fpl,i):
	fi:
	return S1:
else
	ap1 := apsq(fpl,sq[1]):
	ap2 := apsq(fpl,sq[2]):
	fs1 := [seq(newaltfpl(fpl,ap1[j]),j=1..nops(ap1))]:
	fs2 := [seq(newaltfpl(fpl,ap2[j]),j=1..nops(ap2))]:
	for j from 1 to nops(fs1) do
		if fplconn(fs1[j])=iconn then
			S1 := S1 union {fs1[j]}:
		fi:
	od:
	for j from 1 to nops(fs2) do
		if fplconn(fs2[j])=iconn then
			S2 := S2 union {fs2[j]}:
		fi:
	od:
	if nops(S1) = 1 and nops(S2) = 1 then
		print(`Houston`,fpl,i):
		return S1 union S2:
	elif nops(S1) = 1 then
		return S1:
	elif nops(S2) = 1 then
		return S2:
	elif nops(S1) = 0 and nops(S2) = 0 then
		print(`nothing`,fpl,i):
		return S1 union S2:
	else
		print(`Both`,fpl,i):
		return S1 union S2:
	fi:
fi:


end:




#invfplalt(fpl) lists the inverse of the alt paths which can be got from the given fpl
invfplalt := proc(fpl) local i,j,k,B,W,S1,S2,S,n,ipt,ip1,conn,iconn,ap1,ap2,stpart,sq,fs1,fs2,i2,c2:

n := nops(fpl):
conn := fplconn(fpl):
iconn := invei(conn):

i2 := {seq(iconn[i][1],i=1..nops(iconn))}:
c2 := [seq({},i=1..nops(i2))]:

for i from 1 to nops(iconn) do
	for j from 1 to nops(i2) do
		if iconn[i][1] = i2[j] then
			c2[j] := c2[j] union {iconn[i][2]}:
		fi:
	od:
od:
c2 := [seq(c2[i] minus {conn},i=1..nops(c2))]:

#return i2,c2:

S := {seq(seq({[j,k],[j+1,k]},j=0..n),k=1..n)}
union {seq(seq({[j,k],[j,k+1]},k=0..n),j=1..n)}:

B := fpltograph(fpl)[2]:
W := S minus B:

S1 := {}:
S2 := {}:

for i from 1 to nops(i2) do
	if i2[i] <= floor(n/2) then
		ipt := [1,2*i2[i]]:
	elif i2[i] <= n then
		ipt := [2*i2[i]-n,n]:
	elif i2[i] <= n + floor(n/2) then
		ipt := [n,3*n+1-2*i2[i]]:
	else
		ipt := [4*n+1-2*i2[i],1]:
	fi:

	if i2[i] < floor(n/2) then
		stpart := [ipt,[1,ipt[2]+1],[1,ipt[2]+2]]:
		sq := {ipt,[1,ipt[2]+1]}:
	elif i2[i] = floor(n/2) and n mod 2 = 1 then
		stpart := [ipt,[1,ipt[2]+1],[2,ipt[2]+1]]:
		sq := {ipt}:
	elif i2[i] = floor(n/2) and n mod 2 = 0 then
		stpart := [[1,ipt[2]-1],ipt,[2,ipt[2]]]:
		sq := {[1,ipt[2]-1]}:
		ip1 := [1,ipt[2]-1]:
	elif i2[i] < n then
		stpart := [ipt,[ipt[1]+1,n],[ipt[1]+2,n]]:
		sq := {[ipt[1],n-1],[ipt[1]+1,n-1]}:
	elif i2[i] = n  then
		stpart := [[n-1,n],ipt,[n,n-1]]:
		sq := {[n-1,n-1]}:
		ip1 := [n-1,n]:
	elif i2[i] < n+floor(n/2) then
		stpart := [ipt,[n,ipt[2]-1],[n,ipt[2]-2]]:
		sq := {[n-1,ipt[2]-1],[n-1,ipt[2]-2]}:
	elif i2[i] = n+floor(n/2) and n mod 2 = 0 then
		stpart := [[n,2],ipt,[n-1,1]]:
		sq := {[n-1,1]}:
		ip1 := [n,2]:
	elif i2[i] = n+floor(n/2) and n mod 2 = 1 then
		stpart := [ipt,[n,1],[n-1,1]]:
		sq := {[n-1,1]}:
	elif i2[i] < 2*n then
		stpart := [ipt,[ipt[1]-1,1],[ipt[1]-2,1]]:
		sq := {[ipt[1]-1,1],[ipt[1]-2,1]}:
	elif i2[i] = 2*n  then
		stpart := [[2,1],ipt,[1,2]]:
		sq := {ipt}:
		ip1 := [2,1]:
	fi:


	ap1 := {seq(op(apsq(fpl,sq[j])),j=1..nops(sq))}:
	fs1 := {seq(newaltfpl(fpl,ap1[j]),j=1..nops(ap1))}:
	for j from 1 to nops(fs1) do
		if evalb(fplconn(fs1[j]) in c2[i]) then
			#S1 := S1 union {fs1[j]}:
			S1 := S1 union {[i2[i],fs1[j]]}:
		fi:
	od:

od:

return S1:

end: