# Matthew Russell
# Experimental Math
# Homework 20
# I give permission to post this

with(combinat):

read `public_html/em12/C20.txt`:
read `public_html/em12/hw18.txt`:

#####################################################################################

# It is clear that every word of length ≥ k in the alphabet {0,1} must contain one of
# the 2^k patterns of length k (forget about Blanks in this problem). Call a set of
# length-k patterns unavoidable if every sufficiently long word in {0,1} must contain
# one of its members. By using GFunT(A,B,FP,T), with A={0,1}, and FP ranging over
# all 2^(2^k) possible sets of length-k B-free patterns, and see which of them outputs
# a polynomial in T (i.e. a rational function in T whose denominator is 1), write a program
# UA(k), that inputs a positive integer k, and outputs the set of sets of patterns as described.
# Of course for k > 4 is would be hopeless. For example
# UA(1);
# should return:
# {{[0],[1]}}

UA:=proc(k) local A, i, S, R, s:

A:=AllWords({0,1},k):
S:=powerset(A):
R:={}:
for s in S do
	if degree(denom(normal(GFunT({0,1},B,s,T))),T)=0 then
		R:=R union {s}:
	fi:
od:

return R:

end:

# Note: seq(nops(UA(k)),k=1..4) is 1, 3, 39, 8023 - not in Sloane.

# I originally interpreted this to look at all possible sets of B-free patterns
# of length k or shorter. Here is code to do that problem (interesting in its own right).

UAdiff:=proc(k) local A, i, S, R, s:

A:={}:
for i from 1 to k do
	A:=A union AllWords({0,1},i):
od:
S:=powerset(A):
R:={}:
for s in S do
	if degree(denom(normal(GFunT({0,1},B,s,T))),T)=0 then
		R:=R union {s}:
	fi:
od:

return R:

end:


#####################################################################################

# Call a set of patterns of length k almost unavoidable if asymptotically, the number
# of words of length n avoiding that set of patterns has polynomial growth (as opposed
# to exponetial, growth). By using GFunT(A,B,FP,T), with A={0,1}, and FP ranging over
# all 2^(2^k) possible sets of B-free length-k patterns, and see which of them outputs
# a rational function whose denominator has roots whose absolute value is ≥ 1 write a
# program AUA(k), that inputs a positive integer k, and outputs the set of sets of patterns
# as described. Of course for k > 4 is would be hopeless. For example
# AUA(1);
# should return:
# { {[0],[1]}, {[0]}, {[1]} }.

AUA:=proc(k) local A, i, S, R, s, sols, sol, flag:

A:=AllWords({0,1},k):
S:=powerset(A):
R:={}:
for s in S do
	sols:={solve(denom(normal(GFunT({0,1},B,s,T))),T)}:
	flag:=true:
	for sol in sols while flag do
		if abs(evalf(evalc(sol)))<1 then
			flag:=false:
		fi:
	od:
	if flag then
		R:=R union {s}:
	fi:
od:

return R:

end:

# Note: seq(nops(AUA(k)),k=1..4) is 3, 13, 205,  - not in Sloane.

#####################################################################################

# Write a procedure
# IsSubPattern(pat1,pat2,B),
# that inputs two patterns, pat1, pat2 and outputs true if and only if any occurrence
# of pat1 implies an occurence of pat2.
# For example, IsSubPattern([1,1,1],[1,1],B);
# and
# IsSubPattern([1,3,1],[1,B,1],B);
# should both return "true", while
# IsSubPattern([2,1,1],[1,B,1],B);
# should return
# False

IsSubPattern:=proc(pat1,pat2,B) local n1, n2, flag1, flag2, i, offset:

n1:=nops(pat1):
n2:=nops(pat2):
flag1:=false:
for offset from 0 to n1-n2 while not flag1 do
	flag2:=true:
	for i from 1 to n2 while flag2 do 
		if pat1[i+offset]<>pat2[i] and pat2[i]<>B then
			flag2:=false:
		fi:
	od:
	flag1:=flag2:
od:
return flag1:

end:

#####################################################################################

# Using the above, write a procedure
# MinimalFP(B,FP)
# that inputs a symbol for blank, B, and a set of patterns (expressed as lists), and
# returns a subset of FP, whose avoidance implies the avoidance of all the other patterns.
# For example,
# MinimalFP(B,{[1,1,1],[1,1,1,1],[1,B,1],[1,2,1]});
# should return
# {[1,B,1]}

MinimalFP:=proc(B,FP) local F, Fnew, f, f1, flag:

F:=FP:
Fnew:={}:
for f in F do
	flag:=true:
	for f1 in F union Fnew minus {f} while flag do
		if IsSubPattern(f,f1,B) then
			flag:=false:
		fi:
	od:
	if flag then
		Fnew:=Fnew union {f}:
	fi:
	F:=F minus {f}:
od:

return Fnew:

end:

#####################################################################################

# A 2D self avoiding walk (henceforth saw) is a word (list) in the alphabet {-1,1,-2, 2}
# such that in no non-empty consecutive subword ("factor") there are as many 1's and -1's
# AND as many 2's and -2's. For example
# [1,2,-1,-1,-2]
# is a self-avoiding walks, but
# W:=[-1,1,2,-1,-2,-2]
# is not since [op(2..5,W)]=[1,2,-1,-2] has as many 1's as -1's AND as many 2's as -2's.
# [To visualize a self-avoiding walk, think of 1 as "one unit step East",
# -1 as "one unit step West", 2 as "one unit step North", -2 as "one unit step South",
# and start walking from the origin. The condition is equivalent to all the visited places
# being distinct (i.e. "it avoids itself").
#
# Write a recursive "smart brute-force" program SAW(n) that inputs a non-negative integer n,
# and outputs the set of all self-avoiding walks of length n. For example,
# SAW(2);
# should yield
# {[1,1],[1,2],[1,-2],[2,2],[2,1],[2,-1], [-1,-1],[-1,-2],[-1,2],[-2,-2],[-2,-1],[-2,1]}
#
# [Hint: Do it recursively, and for each member of SAW(n-1) find out which of the four possible extensions are "legal"]

SAW:=proc(n) local S1, S, s:
option remember:

if n=0 then
	return {[]}:
fi:

S1:=SAW(n-1):
S:={seq([op(s),-1],s in S1),seq([op(s),1],s in S1),seq([op(s),-2],s in S1),seq([op(s),2],s in S1)}:
for s in S do
	if SAWendcheck(s) then
		S:=S minus {s}:
	fi:		
od:
return S:

end:

SAWendcheck:=proc(w) local flag, i:
option remember:
flag:=true:
for i from 1 to nops(w)/2 while flag do
	if SAWfactorcheck2([op(nops(w)-2*i+1..nops(w),w)]) then
		flag:=false:
	fi:
od:
return not flag:

end:

SAWfactorcheck2:=proc(f) local f1, negone, posone, negtwo, postwo, t:
option remember:

f1:=convert(f,multiset):
negone:=0:
posone:=0:
negtwo:=0:
postwo:=0:
for t in f1 do
	if t[1]=-1 then
		negone:=t[2]:
	elif t[1]=1 then
		posone:=t[2]:
	elif t[1]=-2 then
		negtwo:=t[2]:
	elif t[1]=2 then
		postwo:=t[2]:
	fi:
od:

return evalb(negone=posone and negtwo=postwo):

end:

#####################################################################################

# Find the first 10 members of the sequence "number of 2D self-avoiding walks of length n". Is it in Sloane?

# 4, 12, 36, 100, 284, 780, 2172, 5916, 16268, 44100
# This is A001411.

#####################################################################################

# A 2D self avoiding polygon (henceforth sap) is a word (list) in the alphabet {-1,1,2,-2}
# such that in no non-empty proper consecutive subword ("factor") there are as
# many 1's and -1's AND as many 2's and -2's, but the total number of 1's equals
# the total number of -1's AND the total number of 2's equals the total number of -2's.
# For example
# [1,1,2,-1,-1,-2]
# is a self-avoiding polygon, but
# W:=[1,1,1,2,-1,-1,-2,-1]
# is not since [op(2..7,W)]=[1,1,2,-1,-1,-2] has as many 1's as -1's AND as many 2's as -2's.
#
# Write a program SAP(n) that inputs a non-negative integer n, and outputs the set
# of all self-avoiding polygons of length 2n. (of course there are no saps of odd
# lengths (why?)). For example,
# SAP(1);
# should yield
# {[1,-1],[-1,1],[2,-2],[-2,2]}

SAP:=proc(n) local S1, Good, s, negone, posone, negtwo, postwo, t, S:

if n=0 then
	return {[]}:
fi:
if n=1 then
	return {[1,-1],[-1,1],[2,-2],[-2,2]}:
fi:

S1:=SAW(2*n-1):
Good:={}:
for s in S1 do
	S:=convert(s,multiset):
	negone:=0:
	posone:=0:
	negtwo:=0:
	postwo:=0:
	for t in S do
		if t[1]=-1 then
			negone:=t[2]:
		elif t[1]=1 then
			posone:=t[2]:
		elif t[1]=-2 then
			negtwo:=t[2]:
		elif t[1]=2 then
			postwo:=t[2]:
		fi:
	od:
	if negone=posone and negtwo-postwo=1 and s[-1]<>-2 then
		Good:=Good union {[op(s),2]}:
	elif negone=posone and negtwo-postwo=-1 and s[-1]<>2 then
		Good:=Good union {[op(s),-2]}:
	elif negone-posone=1 and negtwo=postwo and s[-1]<>-1 then
		Good:=Good union {[op(s),1]}:
	elif negone-posone=-1 and negtwo=postwo and s[-1]<>1 then
		Good:=Good union {[op(s),-1]}:
	fi:
od:

return Good:

end:

# seq(nops(SAP(i)),i=1..5): 4, 8, 24, 112, 560 (not in Sloane, and neither is
# 1, 2, 6, 28, 140: the sequence divided by 4)

#####################################################################################

# A finite-memory with memory 2d is a walk as above where each no consecutive subwalk
# of length 2r is a sap, for r=1,2, ..., d. By using
# GFunT(A,B,FP,T);
# with
# A={1,-1,2,-2}
# and FP the union of SAP(1), SAP(2), ..., SAP(d), write a procedure
# GFsawd(d,T)
# that inputs a pos. integer d, and a variable T, and outputs the rational function in T
# whose Maclaurin coefficient of Tn would give the number of memory-d-self-avoiding walks
# of length n. For example,
# GFsawd(1,T);
# should yield (1+T)/(1-3*T).

GFsawd:=proc(d,T) local A, FP, i:

A:={1,-1,2,-2}:
FP:={}:
for i from 1 to d do
	FP:=FP union SAP(i):
od:
return GFunT(A,B,FP,T):

end:

#####################################################################################

# Find GFsawd(d,T), for d=1,d=2,d=3, and if possible d=4. Find the first 30 terms
# in the taylor expansion, and see whether these sequences are in Sloane. 

# GFsawd(1,T): -(T+1)/(3*T-1).
# Taylor expansion coeffs: 1, 4, 12, 36, 108, 324, 972, 2916, 8748, 26244,
#			   78732, 236196, 708588, 2125764, 6377292, 19131876,
#			   57395628, 172186884, 516560652, 1549681956, 4649045868,
#			   13947137604, 41841412812, 125524238436, 376572715308,
#			   1129718145924, 3389154437772, 10167463313316, 30502389939948,
#			   91507169819844, 274521509459532
# (this is in Sloane - A003946 - associated with this generating function.)

# GFsawd(2,T): -(2*T^2+3*T^3+2*T+1)/(T^3+2*T^2+2*T-1).
# Taylor expansion coeffs: 1, 4, 12, 36, 100, 284, 804, 2276, 6444, 18244,
#			   51652, 146236, 414020, 1172164, 3318604, 9395556,
#			   26600484, 75310684, 213217892, 603657636, 1709061740,
#			   4838656644, 13699094404, 38784563836, 109805973124,
#			   310880168324, 880156846732, 2491880003236, 7054953868260,
#			   19973824589724, 56549436919204
# (not in Sloane)

# GFsawd(3,T): -(8*T^8+3*T^7-T^6+8*T^5+3*T^4+4*T^3+2*T^2+2*T+1)/((T+1)*(T^6+T^3-T^2+3*T-1)).
# Taylor expansion coeffs: 1, 4, 12, 36, 100, 284, 780, 2172, 6028, 16732,
#			   46436, 128892, 357748, 992964, 2756060, 7649700,
#			   21232436, 58932564, 163572700, 454010940, 1260148740,
#			   3497657684, 9708067684, 26945626676, 74790042724,
#			   207586580124, 576175473060, 1599227539468, 4438801793148,
#			   12320298939716, 34196112608188
# (not in Sloane)

# GFsawd(4,T):  -(74*T^21+1+11*T^42+2*T+24*T^43-33*T^41-153*T^17-80*T^40+49*T^34+101*T^36
#	  	+52*T^37-16*T^38-72*T^39-61*T^18-117*T^16+141*T^22-36*T^14-132*T^15+141*T^35
#		-54*T^33+3*T^2+6*T^3+84*T^20-21*T^12-T^13-117*T^32+14*T^5+41*T^7+9*T^6+6*T^4
#		-18*T^8+5*T^10+22*T^11+50*T^23-82*T^24+24*T^25+16*T^27+61*T^29+2*T^26+47*T^28
#		-117*T^30-94*T^31-4*T^19+32*T^44+24*T^45+8*T^46+28*T^9)/
#		(6*T^21-1+T^42+2*T+T^41+9*T^17+3*T^34-T^36-4*T^37+T^18-7*T^16-13*T^22-8*T^14
#		-8*T^15+3*T^35+2*T^33+T^2+2*T^3+4*T^20+T^12-7*T^13+5*T^32+2*T^5-5*T^7-T^6+2*T^4
#		+2*T^8+3*T^10-2*T^11+2*T^23-2*T^24+4*T^25+4*T^27-5*T^29+2*T^26-3*T^28-7*T^30
#		-2*T^31+12*T^19-4*T^9).
# Taylor expansion coeffs: 1, 4, 12, 36, 100, 284, 780, 2172, 5916, 16268,
#			   44660, 122596, 336428, 923316, 2533924, 6954340,
#			   19085796, 52380228, 143755252, 394530428, 1082772220,
#		 	   2971623276, 8155495676, 22382417316, 61427608676,
#			   168585505252, 462675873540, 1269794600060, 3484898215012,
#			   9564157518364, 26248430624188
# (not in Sloane)

#####################################################################################