#1234Avoid.txt
read `Cfinite`:

Help:=proc()
if args=NULL then
	print(`Count1234PlusAvoiders(r,k,n)`):
	print(`Guess1234PlusGeneratingFunction(r,k,t)`):
elif args[1] = 'Count1234PlusAvoiders' then
	print(`Count1234PlusAvoiders(r,k,n) outputs the number of words consisting of letters`):
	print(`1,2,...,n, each used r times which simultaneously avoid the patterns 1234 and 1k(k-1)...2`):
	print(`For Example: Count1234PlusAvoiders(2,5,5);`):
elif args[1] = 'Guess1234PlusGeneratingFunction' then
	print(`Guess1234PlusGeneratingFunction(r,k,t) tries to find a generating function for`):
	print(`the sequence given by Count1234PlusAvoiders(r,k,n) in the variable t.`):
	print(`Gives up if 256 terms of the sequence aren't enough to figure out the function`):
	print(`For Example: Guess1234PlusGeneratingFunction(1,4,t);`):
fi:
end:
	
Remove:=proc(L,elt) local i, output, hold:
hold:=[]:
for i from 1 to nops(L) do
	if L[i]<>elt then
		hold:=[op(hold),L[i]]:
	fi:
od:
hold:
end:

Reduce:=proc(L,red) local output, l, elt:
output:=[]:
for l in L do
	output:=[op(output),[]]:
	for elt in l do
		output[-1] := [op(output[-1]),elt - red]:
	od:
od:
output:
end:

ReduceOne:=proc(L,ind) local Lp:
	Lp:=L:
	Lp[ind] := Lp[ind]-1:
	return(Lp):
end:

#b is defined as the largest index which is not greater than any previous letter
Count1234_132Avoiders1:=proc(r,a,b,L) local Lp,bp,ap,output,i,j:
option remember:
if a = 0 and {op(L)} = {0} then
	return(1):
fi:
if member(-1,L) then
	return(0):
fi:

output:=0:
for i from 1 to min(3,b+1) do
	Lp := L:
	ap := a:
	Lp[i] := Lp[i] - 1:
	bp := min(b,i):
	if Lp[i] = 0 then
		for j from i-1 to 1 by -1 do
			Lp[j+1] := Lp[j]:
		od:
		if a > 0 then
			Lp[1] := r:
			ap := ap - 1:
		else
			Lp[1] := 0:
		fi:
		if b < i then
			bp := b+1:
		fi:
	fi:
	output:=output + Count1234_132Avoiders1(r,ap,bp,Lp)
od:
#print(r,a,b,L,output):
output:
end:

Count1234_132Avoiders:=proc(r,n)
	if n < 3 then
		return(mul(binomial(n*r-i*r,r),i=0..(n-1))):
	fi:
	return(Count1234_132Avoiders1(r,n-3,3,[r,r,r])):
end:

NonzeroInds:=proc(L) local i, output:
output:=[]:
for i from 1 to nops(L) do
	if L[i] <> 0 then
		output:=[op(output),i]:
	fi:
od:
output:
end:


Check321:=proc(L) local biggest, bigger,i:
biggest:=0:
bigger:=0:
for i from 1 to nops(L) do
	if L[i] >= biggest then
		biggest := L[i]:
	elif L[i] >= bigger then
		bigger := L[i]:
	else
		return(true):
	fi:
od:
return(false):
end:

Check_k21:=proc(L,k) local biggest,i,j:
biggest := [0$(k-1)]:
for i from 1 to nops(L) do
	for j from nops(biggest) to 1 by -1 do
		if L[i] >= biggest[j] then
			biggest[j] := L[i]:
			break:
		fi:
		if j = 1 then
			return(true):
		fi:
	od:
od:
return(false):
end:

my_min:=proc(li)
	if nops(li) = 0 then
		return(6):
	fi:
	return(min(li)):
end:

RemoveOldInd := proc(r,a,Acted,actors,L,oldind) local ap,Actedp,actorsp,Lp,i,l:
	Lp := L:
	ap := a:
	Actedp := Acted:
	actorsp := actors:
	for i from oldind to 2 by -1 do
		Lp[i] := Lp[i-1]:
	od:
	if ap > 0 then
		Lp[1] := r:
		ap := ap -1:
	else
		Lp[1] := 0:
	fi:
	for i from 1 to nops(Actedp) do
		for l from 1 to nops(Actedp[i]) do
			if Actedp[i][l] < oldind then
				Actedp[i][l] := Actedp[i][l] + 1:
			fi:
		od:
	od:
	for i from 1 to nops(actorsp) do
		if actorsp[i] < oldind then
			actorsp[i] := actorsp[i] + 1:
		fi:
	od:
	return(ap,Actedp,actorsp,Lp):
end:

Count1234_1432Avoiders:=proc(r,n) local max_total:
	max_total:=14:
	if n >= max_total then
		return(Count1234_1432Avoiders1(r,n-max_total,[],[],[r$max_total])):
	else
		return(Count1234_1432Avoiders1(r,0,[],[],[0$(max_total-n),r$n])):
	fi:
end:

Count1234_1432Avoiders1:=proc(r,a,Acted,actors,L) local ap,output,i,j,l,Actedp,actorsp,Lp,additions,viable,deletable_count,available,oldactor,fixed,max_total:
option remember:

max_total := 14:
if a = 0 and {op(L)} = {0} then
	return(1):
fi:
if member(-1,L) then
	return(0):
fi:
output:=0:
available := {op(NonzeroInds(L))}:
if nops(available) > 5 then
	available := {op(-5..-1,available)}:
fi:
for i in available while i <= min(actors) do
	Lp:=L:
	ap:=a:
	Actedp := Acted:
	actorsp := actors:
	Lp[i] := Lp[i] -1:
	if i < min(actors) then
		if nops(actors) > 0 and Acted[-1] = [] then
			actorsp := [op(1..-2, actorsp)]:
			Actedp := [op(1..-2, Actedp)]:
			actorsp := [op(actorsp), i]:
			Actedp := [op(Actedp), []]:
		else
			actorsp := [op(actorsp), i]:
			Actedp := [op(Actedp), []]:
		fi:
	fi:
	for j from 1 to max_total do
		if Lp[j] = 0 and not member(j,{seq(op(Actedp[l]),l=1..nops(Actedp))} union {op(actorsp)}) then
			ap,Actedp,actorsp,Lp := RemoveOldInd(r,ap,Actedp,actorsp,Lp,j):
		fi:
	od:
	output := output + Count1234_1432Avoiders1(r,ap,Actedp,actorsp,Lp):
od:

additions := {}:
for i in available do
	if i > min(actors) and ( i <= min(Acted[1]) or i = max(available)) then
		additions := additions union {i}:
	fi:
od:

for i in additions do
	Lp := L:
	ap := a:
	Actedp := Acted:
	actorsp := actors:
	viable := true:
	deletable_count := 0:
	for j from 1 to nops(Actedp) do
		if i <= actors[j] then
			fixed := sort(available intersect {seq(l, l=actorsp[j]+1..max(available))}):
			if Check321([op(Actedp[j]), seq(fixed[-l], l=1..nops(fixed))]) then
				viable := false:
			fi:
			if nops(available intersect {seq(l, l=actorsp[j]+1..max(available))}) >= 3 then
				viable := false:
			fi:
			deletable_count := deletable_count + 1:
		else
			if Check321([op(Actedp[j]),i]) then
				viable := false:
			fi:
		fi:
	od:
	if viable then
		actorsp:= [op(1+deletable_count..-1, actorsp)]:
		Actedp := [op(1+deletable_count..-1, Actedp)]:
		for j from 1 to nops(Actedp) do
			Actedp[j] := [op(Actedp[j]), i]:
		od:
		Lp[i] := Lp[i] -1:
		for j from 1 to max_total do
			if Lp[j] = 0 and not member(j,{seq(op(Actedp[l]),l=1..nops(Actedp))} union {op(actorsp)}) then
				ap,Actedp,actorsp,Lp := RemoveOldInd(r,ap,Actedp,actorsp,Lp,j):
			fi:
		od:
		output := output + Count1234_1432Avoiders1(r,ap,Actedp,actorsp,Lp):
	fi:
od:
return(output):
end:


Count1234PlusAvoiders:=proc(r,k,n) local num_total,num_available:
	num_available := (k-2)*(2)+1:
	num_total := 2*(num_available-1)+1+num_available:
	if n >= num_total then
		return(Count1234PlusAvoiders1(r,k,n-num_total,[],[],[r$num_total])):
	else
		return(Count1234PlusAvoiders1(r,k,0,[],[],[0$(num_total-n),r$n])):
	fi:
end:

Count1234PlusAvoiders1:=proc(r,k,a,Acted,actors,L) local ap,output,i,j,l,Actedp,actorsp,Lp,additions,viable,deletable_count,available,fixed,num_available,num_total:
option remember:

#if iters =1 then
#	print(r,k,a,Acted,actors,L):
#fi:

num_available := (k-2)*(2)+1:
num_total := 2*(num_available-1)+1+num_available:

if a = 0 and {op(L)} = {0} then
	return(1):
fi:
if member(-1,L) then
	return(0):
fi:
output:=0:
available := {op(NonzeroInds(L))}:
if nops(available) > num_available then
	available := {op(-num_available..-1,available)}:
fi:
for i in available while i <= min(actors) do
	Lp:=L:
	ap:=a:
	Actedp := Acted:
	actorsp := actors:
	Lp[i] := Lp[i] -1:
	if i < min(actors) then
		if nops(actors) > 0 and Acted[-1] = [] then
			actorsp := [op(1..-2, actorsp)]:
			Actedp := [op(1..-2, Actedp)]:
			actorsp := [op(actorsp), i]:
			Actedp := [op(Actedp), []]:
		else
			actorsp := [op(actorsp), i]:
			Actedp := [op(Actedp), []]:
		fi:
	fi:
	for j from 1 to num_total do
		if Lp[j] = 0 and not member(j,{seq(op(Actedp[l]),l=1..nops(Actedp))} union {op(actorsp)}) then
			ap,Actedp,actorsp,Lp := RemoveOldInd(r,ap,Actedp,actorsp,Lp,j):
		fi:
	od:
	output := output + Count1234PlusAvoiders1(r,k,ap,Actedp,actorsp,Lp):
od:

additions := {}:
for i in available do
	if i > min(actors) and ( i <= min(Acted[1]) or i = max(available)) then
		additions := additions union {i}:
	fi:
od:

for i in additions do
	Lp := L:
	ap := a:
	Actedp := Acted:
	actorsp := actors:
	viable := true:
	deletable_count := 0:
	for j from 1 to nops(Actedp) do
		if i <= actors[j] then
			fixed := sort(available intersect {seq(l, l=actorsp[j]+1..max(available))}):
			if Check_k21([op(Actedp[j]), seq(fixed[-l], l=1..nops(fixed))],k-1) then
				viable := false:
			fi:
			if nops(available intersect {seq(l, l=actorsp[j]+1..max(available))}) >= k-1 then
				viable := false:
			fi:
			deletable_count := deletable_count + 1:
		else
			if Check_k21([op(Actedp[j]),i],k-1) then
				viable := false:
			fi:
		fi:
	od:
	if viable then
		actorsp:= [op(1+deletable_count..-1, actorsp)]:
		Actedp := [op(1+deletable_count..-1, Actedp)]:
		for j from 1 to nops(Actedp) do
			Actedp[j] := [op(Actedp[j]), i]:
		od:
		Lp[i] := Lp[i] -1:
		for j from 1 to num_total do
			if Lp[j] = 0 and not member(j,{seq(op(Actedp[l]),l=1..nops(Actedp))} union {op(actorsp)}) then
				ap,Actedp,actorsp,Lp := RemoveOldInd(r,ap,Actedp,actorsp,Lp,j):
			fi:
		od:
		output := output + Count1234PlusAvoiders1(r,k,ap,Actedp,actorsp,Lp):
	fi:
od:
return(output):
end:
	
			

#Guess1234PlusGeneratingFunction(r,k,t) tries to find a generating function for the sequence given by Count1234PlusAvoiders(r,k,n) in the variable t. Gives up if 256 terms of the sequence aren't enough to figure out the function
Guess1234PlusGeneratingFunction:=proc(r,k,t) local num_terms, se, re:
num_terms:=8:
while num_terms <= 512 do
	se:=[seq(Count1234PlusAvoiders(r,k,n),n=0..num_terms)]:
	re:=GuessRec(se):
	if re <> FAIL then
		break:
	fi:
	num_terms:=2*num_terms:
od:
if re = FAIL then
	return(FAIL):
fi:
CtoR(re,t):
end:

Guess1234_132GeneratingFunction:=proc(r,t) local num_terms, se, re:
num_terms:=32:
while num_terms <= 512 do
	se:=[seq(Count1234_132Avoiders(r,n), n=0..num_terms)]:
	re:=GuessRec(se):
	if re <> FAIL then
		break:
	fi:
	num_terms:=2*num_terms:
od:
if re = FAIL then
	return(FAIL):
fi:
CtoR(re,t):
end: