with(combinat, permute):

print(`Type in Help() for help`);

PRINT_RESULT := true:

#####################################################################################
# Print msg if PRINT_RESULT is true or printlevel > 1
#
# msg is a string
#####################################################################################

PrintResult := proc()
	local msg, i;
	
	if PRINT_RESULT and printlevel > 1 then
		msg := ``;
		
		for i from 1 to nargs do
			msg := cat(msg, args[i]);
		od:
		
		print(msg);
	fi;
end:

###################################################################################
# given a list list1, find out if the list is a sublist of list2, starting from start
# return the index of the first letter of the first occurence, 0 otherwise
###################################################################################

IsMemberOf := proc(list1, list2, start)
	local len1, len2, index1, index2;
	
	len1 := nops(list1);
	len2 := nops(list2);
	index1 := start;
	
	if nops(list1) = 0 then 
		RETURN(0);
	fi:
	
	while member(list1[1], list2[index1 .. len2], index2) do
		index1 := index1 + index2 - 1; 

		if index1 + len1 - 1 > len2 then
			RETURN(0);
		fi;

		if list2[index1 .. index1 + len1 - 1] = list1 then
			RETURN(index1);
		fi:
		
		index1 := index1 + 1;
		unassign('index2');
	od;
	
	0;
end:


###################################################################################
# GetBadWords(numofletters, patterns, solidbadwords, n): set of words of length n
# in the alphabet {1, 2, .., numofletters} that have substrings 
# that are of the type of the patterns
# It takes a fourth "hidden" parameter, 
# if it exists and is true, the function will return all goodwords that end with 1, badwords
# otherwise, the function returns only the badwords, 
# which do not include the ones with another badword as a substring
###################################################################################

GetBadWords := proc(numofletters, patterns, solidbadwords, n)
	local bgetgoodwords, ret, formalargs: 

	formalargs := nops([op(1, eval(GetBadWords))]);
	if nargs > formalargs and args[formalargs + 1] = true then
		bgetgoodwords := true;
	else
		bgetgoodwords := false;
	fi:

	ret := RecursivelyGetBadWords(numofletters, patterns, solidbadwords, n, bgetgoodwords);

	if bgetgoodwords then
		PermuteResult(ret[1], numofletters), ret[2];
	else
		ret;
	fi;
end:

RecursivelyGetBadWords := proc(numofletters, patterns, solidbadwords, n, bgetgoodwords)
	local goodwords, badwords, i, j, k, words, word, word1, newpatterns: 
	option remember;

PrintResult(`Calculating result for bad words with length `, n);

	if n = 1 then
		if bgetgoodwords then
#			RETURN([seq([i], i = 1 .. numofletters)], []):
			RETURN([[1]], []):
		else
			RETURN([]);
		fi:

	fi:
	
	newpatterns := OptimizePatterns(patterns);
	
	words := RecursivelyGetBadWords(numofletters, newpatterns, solidbadwords, n - 1, true):
	badwords := words[2];
	words := words[1]; # words are the old good words
	goodwords := []:
	
	if nops(words) = 0 then
		PrintResult(`No good words available for the pattern set: `, newpatterns);

		if bgetgoodwords then
			RETURN([], badwords);
		else
			RETURN(badwords);
		fi:
	else
		# given a set of good words, add one letter at the beginning, and check
		# if the new word is still good
		for i from 1 to nops(words) do
			word := op(i, words):
			
			# for each good word, add an extra letter to it,
			# then check if it is still good
			# we do not have to check bad words, because they are bad already
			# and any thing generated from it will be bad too
			for j from 1 to numofletters do
				word1 := [j, op(word)]:	
				
				# if part of the word is a bad word, ignore it
				for k from 1 to nops(badwords) do
					# only check the head of the word, 
					# because of the way we generated the new word
					if word1[1..nops(badwords[k])] = badwords[k] then
						break;
					fi:
				od:
				if k <= nops(badwords) then
					next;
				fi:

				# determine if it is a good word or bad word
				if IsGoodWordFromStart(word1, newpatterns, solidbadwords, true) then
					if bgetgoodwords then
						goodwords := [op(goodwords), word1]:
					fi:
				else
					badwords := [op(badwords), word1];
				fi:
			od:
		od:
	fi:
	
	badwords := PermuteResult(badwords, numofletters);

PrintResult(`Calculated result for bad words with length `, n);
	
	if bgetgoodwords then
		goodwords, badwords:
	else
		badwords;
	fi:
end:

###################################################################################
# Given a set of patterns, make sure none of them is of the kind of another
# returns the optimized set of pattern to be used by other functions
# The purpose of the function is to provide better performance for the program
###################################################################################

OptimizePatterns := proc(patterns)
	local i, j, newpatterns, replacelist, len, ii;
	
	newpatterns := patterns;
	replacelist := [seq(i, i = 1 .. max(1, seq(op(patterns[ii]), ii = 1 .. nops(patterns))))];
	
	i := 1;
	j := 1;
	len := nops(newpatterns);
	
	while i <= len do
		while j <= len do
			if j <> i then
				if IsWordOfPattern(newpatterns[i], newpatterns[j], replacelist, true) then
					break;
				fi;
			fi:

			j := j + 1;
		od;
		
		# remove any pattern that is of type of another one
		if j <= len then
			newpatterns := subsop(i = NULL, newpatterns);
			len := len - 1;
			i := i - 1;
		fi;
		
		j := 1;
		i := i + 1;
	od;
	
	newpatterns;
end:

###################################################################################
# Check if a word is of the types of patterns, i.e. not a good word
# return false if it is, true otherwise
###################################################################################
 
IsGoodWordFromStart := proc(word, patterns, solidbadwords, bignoretail)
	local i, ii, replacelist:

	if WordHasBadpart(word, solidbadwords) then
		RETURN(false);
	fi;
						
	replacelist := [seq(i, i = 1 .. max(1, seq(op(patterns[ii]), ii = 1 .. nops(patterns))))];
	
	for i from 1 to nops(patterns) do
		if IsWordOfPattern(word, patterns[i], replacelist, bignoretail) then
		    RETURN(false):
		fi:
	od:

	true:
end:

###################################################################################
# Given a pattern, check rcursively if word is of type pattern
# return true if it is, false otherwise
#
# replacelist is the known replacement for symbols for the 
# recursive function, providing the item in the corresponding 
# position is a list
#
###################################################################################

IsWordOfPattern := proc(word, pattern, replacelist, bignoretail)
	local i, head, newpattern, newpattern1, replacelist1, newword, replacement:

	newpattern := pattern;
	newword := word;

	if newword = [] then
		# if the word reaches the end, the pattern must too	
		if newpattern = [] then
			RETURN(true);
		else
			RETURN(false);
		fi:
	elif newpattern = [] then
		# if the pattern reaches the end, the word must too, 
		# unless we can ignore extra tail
		if bignoretail then
			RETURN(true):
		else
			RETURN(false);
		fi;
	fi:
	
	# remove at know patterns at the beginning
	while newpattern <> [] do
		replacement := replacelist[newpattern[1]];
		
		# if the item in the replacement is a list
		# and the list coincide with the beginning 
		# of the word, make the replacement, and 
		# continue recursively.
		# Otherwise, return false.
		if type(replacement, list) then
			if nops(newword) >= nops(replacement) then
				if newword[1..nops(replacement)] = replacement then
					newpattern := subsop(1=NULL, newpattern);
					newword := newword[nops(replacement) + 1..nops(newword)];
					
					if newword = [] then
						if newpattern = [] then
							RETURN(true);
						else
							RETURN(false);
						fi:
					elif newpattern = [] then
						if bignoretail then
							RETURN(true):
						else
							RETURN(false);
						fi;
					fi:
				else
					RETURN(false);
				fi:
			else
				RETURN(false);
			fi:
		else
			break;
		fi:
	od:
	
	# guess the first item in the remaining pattern
	for i from 1 to GetMaxLengthOfNextItem(newword, newpattern, replacelist) do
		head := [op(1 .. i, newword)]:
		
		replacelist1 := subsop(newpattern[1] = head, replacelist);
		
		if IsWordOfPattern(newword, newpattern, replacelist1, bignoretail) then
			RETURN(true);
		fi:
	od:
	
	false;
end:

###################################################################################
# get the maximun possible length of the next item in pattern list within the word
# assuming the item is not available in replacelist yet.
# The purpose of the function is to provide better performance for the program
###################################################################################

GetMaxLengthOfNextItem := proc(word, pattern, replacelist)
	local i, j, len, replacement, num;
	
	len := nops(word);
	num := 0;
	
	if pattern = [] then 
		RETURN(0);
	fi:
	
	# remove the length of all replaced items first
	for i from 1 to nops(pattern) do
		replacement := replacelist[pattern[i]];
		
		if type(replacement, list) then
			len := len - nops(replacement);
		elif pattern[i] <> pattern[1] then
			len := len - 1;
		else
			num := num + 1;
		fi:
	od:
	
	# devide the remaining number by the number of time the item appears in the pattern
	iquo(len, num);
end:

###################################################################################
# Given a list of strings, return the result under all the permutations of numofletters
###################################################################################

PermuteResult := proc(resultlist, numofletters)
	local i, j, k, len, ret, perm, lenperm, list1, list2, len2 ;

	len := nops(resultlist);
	perm := permute(numofletters);
	lenperm := nops(perm);
	ret := {};
	
	for i from 1 to len do
		list1 := resultlist[i];
		len2 := nops(list1);

		for j from 1 to lenperm do
			list2 := [];
		
			for k from 1 to len2 do
				list2 := [op(list2), perm[j][list1[k]]];
			od;
			ret := ret union {list2};
		od;
	od;
	
	convert(ret, list);
end:

###################################################################################
# Given a list of strings, return all the strings that do not contain badstrings,
# which are solid words
###################################################################################

RemoveBadWords := proc(wordlist, badwords)
	local i, j, lenlist, lenwords, ret;
	
	ret := [];
	lenlist := nops(wordlist);
	
	for i from 1 to lenlist do
		if not WordHasBadpart(wordlist[i], badwords) then
			ret := [op(ret), wordlist[i]];
		fi;
	od;
	
	ret;
end:

###################################################################################
# Given a string, check if it contains certain bad strings,
# which are solid words
###################################################################################

WordHasBadpart := proc(word, badwords)
	local j, lenwords;

	lenwords := nops(badwords);

	for j from 1 to lenwords do
		if IsMemberOf(badwords[j], word, 1) <> 0 then
			RETURN(true);
		fi;
	od;
	
	false;
end:


###################################################################################
# Versify a homomorphism is good by determine the images under the homomorphism
# do not has instances of the patterns
#
# startword is the word to be started with for the conversion
# numofletters is the number of letters in the alphabet
# conversion is the homomorphism 
# patterns is a list of patterns to be checked against for words not to comply with
# N is the minimum length of the words to be checked
###################################################################################

IsGoodConversion := proc(startword, numofletters, conversion, patterns, solidbadwords, N)
	local i, j, len, convset;
	
	convset := [op(ConvertString(startword, conversion))];
	len := nops(convset);
	
	for i from 1 to len do
		if nops(convset[i]) < N then
			if not IsGoodConversion(convset[i], numofletters, conversion, patterns, solidbadwords, N) then
				RETURN(false);
			fi;
		else
			if IsBadWord(convset[i], patterns, solidbadwords) then
				RETURN(false);
			fi;
		fi;
	od;
	
	true;
end:

###################################################################################
# string is a list of alphabet
# conversion is of format 
# [ [[input1], [[output11], [output12], ...]], [[input2] , [[output21], [output22], ...]], ...]
###################################################################################

ConvertString := proc(string, conversion)
	local i, j, k, retset, lenconv, lenstring, leninput, result;

	lenconv := nops(conversion);
	lenstring := nops(string):
	retset := {};

	for i from 1 to lenconv do
		leninput := nops(conversion[i][1]);
		
		if leninput <= lenstring and [op(1..leninput, string)] = conversion[i][1] then
			for j from 1 to nops(conversion[i][2]) do
				if leninput = lenstring then
					retset := retset union {op(conversion[i][2])}:
				else
					result := ConvertString([op(leninput + 1 .. nops(string), string)], conversion);
				
					for k from 1 to nops(result) do
						retset := retset union {[op(conversion[i][2][j]), op(result[k])]};
					od;
				fi;
			od;
		fi;
	od;

	retset;
end:

###################################################################################
# Find out if the word has a substring that follows the patterns
# Since IsWordOfPattern can ignore redundant tails, 
# we only check the substrings that start at the middle of the word 
# to the end of the word
#
# word is a list of letters
# patterns is a list of abstract patterns
###################################################################################

IsBadWord := proc(word, patterns, solidbadwords)
	local len, i, ii, j, replacelist, tempword;

	if WordHasBadpart(word, solidbadwords) then
		RETURN(true);
	fi;
	
	len := nops(word);
	replacelist := [seq(i, i = 1 .. max(1, seq(op(patterns[ii]), ii = 1 .. nops(patterns))))];
	
	for j from 1 to len do
		for i from 1 to nops(patterns) do
			if IsWordOfPattern(word[j..len], patterns[i], replacelist, true) then
			    RETURN(true):
			fi:
		od:
	od;
	
	false:
end:



###################################################################################
#
#
# GUESSING TOOLS
#
#
###################################################################################



###################################################################################
# Given a list of words wordlistfrom, try to convert them into another list of 
# words worldlistto, so that words in largewordlistfrom can be successfully and legally 
# converted into words in largewordlistto applying the same conversion 
###################################################################################

DecodeWords := proc(wordlistfrom, wordlistto, largewordlistfrom, largewordlistto)
	local i, j, len, wordcountfrom, wordcountto, count, countlist;
	
	wordcountfrom := GetListHeadTailCount(wordlistfrom, largewordlistfrom);
	wordcountfrom := sort(wordcountfrom, SortWordCount);
	wordcountto := GetListHeadTailCount(wordlistto, largewordlistto);
	wordcountto := sort(wordcountto, SortWordCount);

	RecursivelyDecodeWords(wordcountfrom, wordcountto, largewordlistfrom, 
		largewordlistto, []);
end:

###################################################################################
# Given a list of words with head/tail counts, try to recursively convert the words
# into another set of words, so that when applying the same conversion to 
# words largewordlistfrom, the result will be in largewordlistto
# decoder has the format of [[[word-to-be-decoded], [word-decoded-to]], [[word-to-be-decoded], [word-decoded-to]], ...]
#
# input one decoder, returns a list of possible extensions of decoder
###################################################################################

RecursivelyDecodeWords := proc(wordcountfrom, wordcountto, largewordlistfrom, 
		largewordlistto, decoder)
	local i, j, k, lenfrom, lento, 
		decodefrom, decodeto, lendecoder, newdecoder, newdecoderlist,
		newfrom, newto, index1, index2, lendecoderlist;
		
	lenfrom := nops(wordcountfrom);
	lento := nops(wordcountto):
	lendecoderlist := nops(decoderlist);
	
	newdecoderlist := [];

	# since this is a recursive function and every word in wordcountfrom will 
	# be eventually converted, we can only try to convert the first word 
	# in the list, and let the recursive part to take care of the rest
	i := 1;	
#	for i from 1 to lenfrom do
		for j from 1 to lento do
			# the number of times that a letter appears as head/tail in 
			# the target of the conversion must be greater than or equal to
			# those of the source
			if wordcountto[j][2] >= wordcountfrom[i][2] and
					wordcountto[j][3] >= wordcountfrom[i][3] then
				decodefrom := wordcountfrom[i][1];
				decodeto := wordcountto[j][1];
				
#				if lendecoderlist = 0 then
#					# if there is nothing in the decoderlist yet, make one
#					
#					newdecoder := [[decodefrom, decodeto]];
#				
#					# test if the new deocder is good
#					if not TestDecoder(largewordlistfrom, largewordlistto, newdecoder) then
#						next;
#					fi;
#
#					newfrom := subsop(i = NULL, wordcountfrom);
#					newto := subsop(j = NULL, wordcountto);
#					newdecoder := RecursivelyDecodeWords(newfrom, newto, largewordlistfrom, 
#							largewordlistto, newdecoder);
#					newdecoderlist := [op(newdecoderlist), op(newdecoder)];
#				else
					# try to grow the decoder list
					newdecoder := [op(decoder), [decodefrom, decodeto]];
				
					# test if the new deocder is good
					if not TestDecoder(largewordlistfrom, largewordlistto, newdecoder) then
						next;
					fi;
				
					# since the new decoder is good for now
					# update the wordfrom and wordto lists by erasing the decoded words
					newfrom := subsop(i = NULL, wordcountfrom);
					newto := subsop(j = NULL, wordcountto);
	
					if nops(newfrom) <> 0 then
						newdecoder := RecursivelyDecodeWords(newfrom, newto, largewordlistfrom, 
								largewordlistto, newdecoder);
										
						if newdecoder <> [] then
							newdecoderlist := [op(newdecoderlist), op(newdecoder)];
						fi;
					else
						newdecoderlist := [op(newdecoderlist), op(newdecoder)];
					fi;
#				fi;
			fi;
		od;
#	od;
	
	newdecoderlist;
end:

###################################################################################
# Given a decoder, find out if it can convert words from largewordlistfrom 
# into largewordlistto. words are tested only if it can be assembled by the words 
# appearing decoder.
#
# return true if all such words can be converted into target word set, 
# even if no conversion actually happenned
# false otherwise
###################################################################################

TestDecoder := proc(largewordlistfrom, largewordlistto, decoder)
	local i, j, decoderwordlist, conversion, lenwordfrom, wordto, ret;
	
	if decoder = [] then 
		RETURN(true);
	fi;
	
	decoderwordlist := [];
	conversion := [];
	
	for i from 1 to nops(decoder) do
		decoderwordlist := [op(decoderwordlist), decoder[i][1]];
		conversion := [ op(conversion), [decoder[i][1], [decoder[i][2]]] ];
	od;
	
	lenwordfrom := nops(decoderwordlist[1]);
	
	for i from 1 to nops(largewordlistfrom) do
		wordto := largewordlistfrom[i];
		
		while wordto <> [] do
			if nops(wordto) >= lenwordfrom then
				if not member(wordto[1..lenwordfrom], decoderwordlist) then
					break;
				else
					wordto := wordto[lenwordfrom + 1 .. nops(wordto)]; 
				fi;
			else
				ERROR(`Unexpected input`);
			fi;
		od;
		
		if wordto <> [] then
			next;
		fi;
		
		ret := ConvertString(largewordlistfrom[i], conversion)[1];
		if not member(ret, largewordlistto) then
			RETURN(false);
		fi;
	od;
	
	true;
end:

###################################################################################
# sort the head/tail counts by head counts first
###################################################################################

SortWordCount := proc(wordcount1, wordcount2)
	if wordcount1[2] > wordcount2[2] then
		true;
	elif wordcount1[2] = wordcount2[2] then
		if wordcount1[3] >= wordcount2[3] then
			true;
		else
			false;
		fi;
	else
		false;
	fi;
end:

###################################################################################
# Given a list, count the times the words in the list appears as prefixes and suffixes 
# in the words in the largewordlist
#
# returns [[word, headcount, tailcount], [word, headcount, tailcount], ...]
###################################################################################

GetListHeadTailCount := proc(wordlist, largewordlist)
	local lenword, lenlist, lentemp, i, j, ret;
	
	lenword := nops(wordlist);
	ret := [];
	
	for i from 1 to lenword do
		ret := [op(ret), [wordlist[i], GetHeadTailCount(wordlist[i], largewordlist)]];
	od;
	
	ret;
end:

###################################################################################
# Given a word, count the times it appears as prefixes and suffixes 
# in the words in the largewordlist
#
# returns [word, headcount, tailcount]
###################################################################################

GetHeadTailCount := proc(word, largewordlist)
	local lenword, lenlist, lentemp, i, j, head, tail;
	
	lenword := nops(word);
	lenlist := nops(largewordlist);
	head := 0;
	tail := 0;
	
	for i from 1 to lenlist do
		lentemp := nops(largewordlist[i]);
		
		if lentemp >= lenword then
			if largewordlist[i][1..lenword] = word then
				head := head + 1;
			fi;
			
			if largewordlist[i][lentemp - lenword + 1..lentemp] = word then
				tail := tail + 1;
			fi;
		fi;
	od;
	
	head, tail;
end:


###################################################################################
#
#
# TESTING TOOLS
#
#
###################################################################################


###################################################################################
# Reverse the order of the letters in a word
###################################################################################

ReverseWord := proc(word)
	local i, len, ret;
	
	len := nops(word);
	ret := [];
	
	for i from 1 to len do
		ret := [word[i], op(ret)];
	od;
	
	ret;
end:

###################################################################################
# Convert all occurance of a letter in a word into another letter
###################################################################################

ConvertLetter := proc(word, letterfrom, letterto)
	local i, len, ret;
	
	len := nops(word);
	ret := word;
	
	for i from 1 to len do
		if ret[i] = letterfrom then
			ret := subsop(i = letterto, ret);
		fi;
	od;
	
	ret;
end:



###################################################################################
#
#
# GUESSING TOOLS
#
#
###################################################################################

###################################################################################
# Main Function
#
# Find all Brinkhuis-pairs of the maxinum possible length 
# when the alphabet has numberofletters letters, and the strings
# are of N letters long and avoiding patterns and solidbadwords.
###################################################################################

GetOptimalTripleSets := proc(numbeofletters, patterns, solidbadwords, N)
	local i, goodwords, goodtriples, triplelist, tripleset;

	goodwords := GetBadWords(numbeofletters, patterns, solidbadwords, N, true)[1]:
	triplelist := GetPermissibleTriples(goodwords, patterns, numbeofletters);
	tripleset := GetInitTripleSets(triplelist):
	GetOptimalTripleSets(tripleset, triplelist);
end:

###################################################################################
# Given a list of triples, find all the triple sets that consist of
# one word in the first two elements and two words in the third
# This will be the base sets to be grown to the optimal triplesets
#
# Assume the triples are well ordered
# #so that the same first elements are always grouped together#
###################################################################################

GetInitTripleSets := proc(triplelist)
	local i, i1, i2, j, j1, j2, k, len, word1, newtriple, newtripleset, newtriplelist, newlen, 
		index1, index2, index3, index, indexlist;
		
	len := nops(triplelist);
	newtripleset := [];
	index := 1;
	
	for i1 from 1 to len do
		word1 := triplelist[i1];
		indexlist := [i1];
		
		for i2 from i1 + 1 to len do
			if member(word1[1], triplelist[i2]) and 
					member(word1[2], triplelist[i2]) then
				indexlist := [op(indexlist), i2];
			else
				break;
			fi;
		od;
		
		if nops(indexlist) > 1 then
			for j1 from 1 to nops(indexlist) do
				for j2 from j1 + 1 to nops(indexlist) do
					newtriple := [[word1[1]], [word1[2]], 
						[triplelist[indexlist[j1]][3], triplelist[indexlist[j2]][3]]];
					newtripleset := [op(newtripleset), 
						[newtriple, 1, [indexlist[j1], indexlist[j2]]]];
				od;
			od;
		end if;
		
		i1 := i2 - 1;
	od;
	
	newtripleset;
end:	

###################################################################################
# Given a list of triplesets and a list triples
# return the set of optimal triplesets that are grown from the input triplesets
#
# Assume the triples are well ordered
# #so that the same first elements are always grouped together#
###################################################################################

GetOptimalTripleSets := proc(triplesetlist, triplelist)
	local i, j, word, newtriple, newtriplesetlist, tripleset, newtripleset,
		bchanged, nfirst, itemlist, lenset, lenlist, ret, nlasti, optimalnumber,
		notavaillist;

	lenset := nops(triplesetlist);
	lenlist := nops(triplelist);
	newtriplesetlist := triplesetlist;
	optimalnumber := 1;
	
	# notavaillist is a list of entries of triplelist 
	# that do not appear in triplesetlist
	# therefore they can be safely ignored
	ret := {}:
	for i from 1 to nops(triplesetlist) do
		ret := ret union convert(triplesetlist[i][3], set);
	od:
	notavaillist := convert({seq(i, i = 1..nops(triplelist))} minus ret, list);
	
	i := 1;
	nlasti := 0;
	while i <= lenset do
		tripleset := newtriplesetlist[i];
		bchanged := false;
		nfirst := 0;

if nlasti - i >= 10 then
PrintResult(`Reached `, i, `th item out of `, nops(newtriplesetlist), 
	` total. Current entry is `, tripleset[3], `. time = `, time());
nlasti := i;
fi;

		for j from tripleset[2] to lenlist do
			if member(j, notavaillist) then
				# if j is in the ignore list, skip it
				next;
			elif member(j, tripleset[3]) then
				# if j is already used by the tripleset, skip it
				next;
			fi;
			
			# test if the current entry in triplelist is compatible with the tripleset
			# return the result set if compatible, otherwise NULL
			ret := MergeTripleIntoTripleSets(triplelist[j], tripleset[1], triplelist);

			if ret <> NULL then
				bchanged := true;
						
				if nfirst = 0 then
					nfirst := j;
					newtripleset := [ret, nfirst + 1, [op(tripleset[3]), j]];
				else
					newtripleset := [ret, nfirst, [op(tripleset[3]), j]];
				fi;
			
				# if the entry list already exists in some other tripleset, maybe not 
				# in the exact order, we can skip the entry, 
				# because the result will be the same
				if DoesHeadExistInTripleSet([op(tripleset[3]), j], newtriplesetlist) then
PrintResult(`Skipping `, i, ` for `, [op(tripleset[3]), j], 
	` because it exists. `, nops(newtriplesetlist));
					next;
				fi;

				newtriplesetlist := [op(newtriplesetlist), newtripleset];

				lenset := lenset + 1;
			fi;
		od;
		
		if bchanged then
			# if the tripleset was used to expand to other triplesets,
			# the current one can be removed
			newtriplesetlist := subsop(i = NULL,  newtriplesetlist);
			
			lenset := lenset - 1;
			i := i - 1;
		else
			# otherwise, this is the last of its kind
			# it shall be removed only when the optimalnumber gets larger
			# because we only want the optimal combination
			ret := min(nops(tripleset[1][1]), nops(tripleset[1][2]), nops(tripleset[1][3]));
			if ret < optimalnumber then
				newtriplesetlist := subsop(i = NULL, newtriplesetlist);
PrintResult(`Tripleset is removed because the optimal number is larger: `, optimalnumber, tripleset[1]);					
				lenset := lenset - 1;
				i := i - 1;
			elif ret > optimalnumber then
				optimalnumber := ret;
PrintResult(`Triplesets are removed because the optimal number grows: `, optimalnumber, `, `, i, `, `, lenset, `, `, nops(newtriplesetlist));

				newtriplesetlist := newtriplesetlist[i..lenset];
				lenset := nops(newtriplesetlist); #lenset - i + 1;
				i := 1;
			elif IsTripleSubset(i, tripleset[3], newtriplesetlist) then
PrintResult(`Tripleset is removed because it is a subset: `, optimalnumber, tripleset[1]);					
				newtriplesetlist := subsop(i = NULL, newtriplesetlist);
					
				lenset := lenset - 1;
				i := i - 1;
			fi;
		fi;
		
		i := i + 1;
	od;
	
	newtriplesetlist;
end:

###################################################################################
# tripleindex is the index of the newly updated triple
# listindex is the list of indexes of triplelist that composite the triple
# tripleset is the list of permissible triples
###################################################################################

DoesHeadExistInTripleSet := proc(listindex, triplesetlist)
	local i, len, tempset, indexset;

	indexset := {op(listindex)};
	len := nops(listindex);
	
	for i from 1 to nops(triplesetlist) do
		if nops(triplesetlist[i][3]) < len then
			next;
		fi;
		
		tempset := {op(triplesetlist[i][3][1..len])};
		
		if tempset = indexset then
			RETURN(true);
		fi;
	od;
	
	false;
end:
	
###################################################################################
# tripleindex is the index of the newly updated triple
# listindex is the list of indexes of triplelist that composite the triple
# tripleset is the list of permissible triples
###################################################################################

IsTripleSubset := proc(tripleindex, listindex, triplesetlist)
	local i, tempset, indexset;
	
	indexset := {op(listindex)};
	
	for i from 1 to tripleindex - 1 do
		tempset := {op(triplesetlist[i][3])};
		
		if indexset minus tempset = {} then
			RETURN(true);
		fi;
	od;
	
	false;
end:

###################################################################################
# Given a triple and a tripleset, merge the two together
# if the result is still a permissible tripleset, return it
# otherwise, return NULL
###################################################################################

MergeTripleIntoTripleSets := proc(triple, tripleset, triplelist)
	local i1, i2, i3, j, newtriple, word;

	for i1 from 1 to nops(tripleset[1]) do
		for i2 from 1 to nops(tripleset[2]) do
			for i3 from 1 to nops(tripleset[3]) do
				word := [triple[1], tripleset[2][i2], tripleset[3][i3]];
				if not member(word, triplelist) then
					RETURN(NULL);
				fi;

				word := [tripleset[1][i1], triple[2], tripleset[3][i3]];
				if not member(word, triplelist) then
					RETURN(NULL);
				fi;

				word := [tripleset[1][i1], tripleset[2][i2], triple[3]];
				if not member(word, triplelist) then
					RETURN(NULL);
				fi;
			od;
		od;
	od;
	
	[ convert(convert(tripleset[1], set) union {triple[1]}, list),
	convert(convert(tripleset[2], set) union {triple[2]}, list), 
	convert(convert(tripleset[3], set) union {triple[3]}, list) ]; 
end:

###################################################################################
# Given a list of words, find out all the triples (A, B, C) in the list that are compatible 
# that is, both the triples are all permissible.
###################################################################################

GetPermissibleTriples := proc(words, patterns)
	local i, j, k, list1, list2, list3, word1, word2, word3, index1, index2, index3, 
		asdf, ret;
	
	ret := [];
	
	list1 := GetAllWordsForFixedPrefix(words, [1, 2], true);
	
	for i from 1 to nops(list1) do
		word1 := list1[i];
		unassign('index1');
		asdf := member(word1, words, index1);
		
		list2 := words[index1 + 1 .. nops(words)];
		list2 := GetPermissibleWords(list2, patterns, list1[i]);
		
		for j from 1 to nops(list2) do
			word2 := list2[j];
			unassign('index2');
			asdf := member(word2, words, index2);
			
			list3 := words[index2 + 1 .. nops(words)];
			list3 := GetPermissibleWords(list3, patterns, list1[i], list2[j]);
			
			for k from 1 to nops(list3) do
				word3 := list3[k];
				
				if IsTriplePermissible(word1, word2, word3, patterns) then	
					ret := [op(ret), [word1, word2, word3]];
				fi;
			od;
		od;
	od;
	
	ret;
end:

###################################################################################
# A triple [u1, u2, u3] is permissible if for any square-free string [w1, w2, w3] 
# of length 3, [u(w1), u(w2), u(w3)] is also square-free
###################################################################################

IsTriplePermissible := proc(word1, word2, word3, patterns)
	local i, j, goodlist, lenlist, conversion, result;

	goodlist := [[1, 2, 3], [2, 1, 3], [2, 3, 1], [1, 3, 2], 
			[1, 2, 1], [3, 1, 2], [3, 2, 1], [1, 3, 1], [3, 2, 3], 
			[3, 1, 3], [2, 3, 2], [2, 1, 2]]:
	lenlist := nops(goodlist);
	conversion := [[[1], [word1]], [[2], [word2]], [[3], [word3]]];
	
	for i from 1 to lenlist do
		result := ConvertString(goodlist[i], conversion)[1];
		if IsBadWord(result, patterns, []) then
			RETURN(false);
		fi;
	od;
	
	true;
end:

###################################################################################
# Given a list of words, find out all words (A) in the list that compatible with
# a list of given words (B), that is, both AB and BA are not in pattern.
# The parameters from the third and on are the list of given words
###################################################################################

GetPermissibleWords := proc(words, patterns, word1)
	local lenwords, lenword1, i, j, k, l, ret, wordindex, 
		word, lenword, argv, lenargv, lenpat, replacelist, ii;
	
	wordindex := nops([op(1, eval(GetPermissibleWords))]) - 1;
	lenwords := nops(words);
	lenword1 := nargs - wordindex;
	ret := [];
	lenpat := nops(patterns);
	
	for i from 1 to lenwords do
		word := words[i];
		lenword := nops(words[i]);
		
		for j from 1 + wordindex to lenword1 + wordindex do
			argv := args[j];
			lenargv := nops(args[j]);
			replacelist := [seq(i, i = 1 .. max(1, seq(op(patterns[ii]), ii = 1 .. nops(patterns))))];

			for k from max(lenargv, lenword) by -1 to 1 do
				for l from 1 to lenpat do
					if k <= lenargv and IsWordOfPattern([op(k..lenargv, argv), op(word)], 
							patterns[l], replacelist, true) then
						break;
					fi;

					if k <= lenword and IsWordOfPattern([op(k..lenword, word), op(argv)], 
							patterns[l], replacelist, true) then
						break;
					fi;
				od;
				
				if l <= lenpat then
					break;
				fi;
			od;

			if k >= 1 then
				break;
			fi;
		od;

		if j > lenword1 + wordindex then
			ret := [op(ret), words[i]];
		fi;
	od;
	
	ret;
end:

###################################################################################
# Given a list of words, get all the words that start with prefix
###################################################################################

GetAllWordsForFixedPrefix := proc(words, prefix, bWithOrWithout)
	local lenwords, lenprefix, i, j, ret;
	
	lenwords := nops(words);
	lenprefix := nops(prefix);
	ret := [];
	
	for i from 1 to lenwords do
		if words[i][1..lenprefix] = prefix then
			if bWithOrWithout then
				ret := [op(ret), words[i]];
			fi;
		else
			if not bWithOrWithout then
				ret := [op(ret), words[i]];
			fi;
		fi;
	od;
	
	ret;
end:

###################################################################################
# Given a list of words, get all the words that start with suffix
###################################################################################

GetAllWordsForFixedSuffix := proc(words, suffix, bWithOrWithout)
	local lenwords, lensuffix, i, j, ret;
	
	lenwords := nops(words);
	lensuffix := nops(suffix);
	ret := [];
	
	for i from 1 to lenwords do
		if words[i][nops(words[i]) - lensuffix + 1..nops(words[i])] = suffix then
			if bWithOrWithout then
				ret := [op(ret), words[i]];
			fi;
		else
			if not bWithOrWithout then
				ret := [op(ret), words[i]];
			fi;
		fi;
	od;
	
	ret;
end:


###################################################################################
# Help function
###################################################################################

Help := proc()
	if args=NULL then
		print(`Contains the following procedures: Main`):
	fi:

	if nops([args])=1 and op(1, [args])=Main then
		print(`Main is to find lowerbound for ternary square-free words`);
		print(`by searching for the Brinkhuis triples. In this program we are looking`);
		print(`for triples whose second and third sets of words are the images of the`);
		print(`permutations (1,2,3) and (1,3,2).`);
		print(``);
		print(`Main function takes one parameter: initwords, which are a list of known`);
		print(`square-free words of the same length a. The function will return`);
		print(`Brinkhuis triples whose words are of length a + 12`);
	fi:
end:

###################################################################################
# Main Function
# 
# Given a set of good words of length n, find the optimal Brinkhuis-pairs
# of length n+6
###################################################################################

Main := proc(initwords)
	local apw, gwft, acp, act, mts;
	
	apw := GetAllPossibleWords(initwords);
	gwft := FindGoodWordsForTriples(apw, [[1, 1]]);
	acp := FindAllCompatiblePairs(gwft);
	act := FindAllCompatibleTriples(gwft, acp);
	mts := GetMaxTripleSet(act, gwft);
end:

###################################################################################
# Get words starting with specified heads and tails, excluding the bad ones
# derived form FindHeadAndTails
###################################################################################

GetAllPossibleWords := proc(initwords)
	local i, j, len, lenword, allwords, lenallwords, 
		lenheadtail, badheadandtail, lenbadhead, lenbadtail;
	
	allwords := [op(ExtendWords(initwords, [1,2,3,1,3,2], [2,3,1,3,2,1], [[1,1]], [])), 
		op(ExtendWords(initwords, [1,2,3,2,1,3], [3,1,2,3,2,1], [[1,1]], []))];
	lenallwords := nops(allwords);
	lenword := nops(allwords[1]);
	
	badheadandtail := FindHeadsAndTails(15)[2];
	lenheadtail := nops(badheadandtail);
	for i from 1 to lenheadtail do
		if nops(badheadandtail[i][1]) > 6 then
			break;
		fi;		
	od;
	badheadandtail := badheadandtail[i .. lenheadtail];
	lenheadtail := nops(badheadandtail);
	
	i := 1;
	while i <= lenallwords do
		for j from 1 to lenheadtail do
			if 	allwords[i][1..nops(badheadandtail[j][1])] = badheadandtail[j][1] and
					allwords[i][lenword - nops(badheadandtail[j][2] + 1 .. lenword)] 
					= badheadandtail[j][2] then
				allwords := subop(i = NULL, allwords);
				i := i - 1;
				lenallwords := lenallwords - 1;
				
				break;
			fi;
		od;
	
		i := i + 1;
	od;
	
	allwords;
end:

###################################################################################
# Given a list of words and fixedwords A and B, find all the words C such that 
# C start with A and end with B
###################################################################################

FindWordsWithFixedStartAndEnd := proc(wordlist, fixedstartword, fixedendword)
	local i, len, reverse, ret, word, lenfixstart, lenfixend, lenword;
	
	len := nops(wordlist);
	lenfixstart := nops(fixedwordstart);
	lenfixend := nops(fixedwordend);
	ret := [];
	
	for i from 1 to len do
		word := wordlist[i];
		lenword := nops(word);
		
		if lenword >= lenfixstart + lenfixend then
			if word[1..lenfixstart] = fixedstartword and 
					word[lenword - lenfixend + 1 .. lenword] = fixedendword then
				ret := [op(ret), word];
			fi;
		fi;
	od;
	
	ret;
end:

###################################################################################
# Given a list of words and a string A and B, find all the words C such that 
# ACB is a good word, and returns ACB
###################################################################################

ExtendWords := proc(wordlist, extensionstart, entensionend, patterns, solidbadwords)
	local i, j, len, temp, ret;
	
	ret := [];
	len := nops(wordlist);
	
	for i from 1 to len do
		temp := [op(extensionstart), op(wordlist[i]), entensionend];
		
		if not IsBadWord(temp, patterns, solidbadwords) then
			ret := [op(ret), temp];
		fi;
	od;
	
	ret;
end:

###################################################################################
# Given a list of words, find all the words A such that 
# [A, A2, A3] and [r(A), r(A2), r(A3)] are comptible permissible triples
# where A2 and A3 are images of A under permutation (2,3,1) and (3,1,2)
# and r(A) is the reverse of A. Comptible means [(A, r(A)), (A2, r(A2), (A3, r(A3))]
# is a Brinkhuis-pair
###################################################################################

FindGoodWordsForTriples := proc(wordlist, patterns)
	local i, j, word, word2, word3, rword, len, ret, ret2, ret3, index;
	
	ret := []:
	len := nops(wordlist);
	
	for i from 1 to len do
		word := wordlist[i];
		word2 := PermuteWord123(word);	
		word3 := PermuteWord132(word);
		rword := ReverseWord(word);
		
		if rword <> word and member(rword, ret) then
			next;
		fi;

		if IsTriplePermissibleForCompatibleWords(word, word2, word3, patterns) then
			# if the word is not its own reverse, then the reverse must be good too
			if rword <> word then
				if IsTriplePermissibleForCompatibleWords(rword, ReverseWord(word2), 
						ReverseWord(word3), patterns) and
						not member(word, ret) and
						AreTwoWordsCompatible(word, rword) then
					ret := [op(ret), word, rword];
				fi;
			else
				ret := [op(ret), word];
			fi;
		fi;
	od;
	
	ret;
end:

###################################################################################
# Test extra cases for 121 and 131, where the first and last words are not the same
# these cases are not checked in the conversions of IsTriplePermissibleForCompatibleWords
###################################################################################

TestExtraCases := proc(word1, word2, word3, patterns)
	local word;
	
	word := [op(word1), op(word2), op(PermuteWord123(word3))];
	if IsBadWord(word, patterns, []) then
		RETURN(false);
	fi;
		
	word := [op(word1), op(word3), op(PermuteWord132(word2))];
	if IsBadWord(word, patterns, []) then
		RETURN(false);
	fi;
	
	true;
end:

###################################################################################
# Verify that triples generated from two words can form a legal Brinkhuis-pair
# Assuming each triple generated from each word is permissible
###################################################################################

AreTwoWordsCompatible := proc(word1, word2)
	local i1, i2, i3, i, len, ret, list1, list2, list3;
	
	list1 := [word1, word2];
	list2 := [PermuteWord123(word1), PermuteWord123(word2)];
	list3 := [PermuteWord132(word1), PermuteWord132(word2)];
	ret :=[]:

	# loop through all possible combinations of triples
	# except the two generated from one word
	for i1 from 1 to 2 do
		for i2 from 1 to 2 do
			for i3 from 1 to 2 do
				if not (i1 = 1 and i2 = 1 and i3 = 1) and
						not (i1 = 2 and i2 = 2 and i3 = 2) and
						not IsTriplePermissibleForCompatibleWords(list1[i1], list2[i2], list3[i3], [[1,1]]) then	
					RETURN(false);
				fi;
			od;
		od;
	od;
	
	true;
end:

###################################################################################
# Given a list of words, find all pairs of indexes of words 
# that are compatible with each other
###################################################################################

FindAllCompatiblePairs := proc(wordlist)
	local ret, ret1, len, i, j;
	
	len := nops(wordlist);
	ret := [];
	
	for i from 1 to len do
print(`FindAllCompatiblePairs: `, i, time()):
		ret1 := [];
		
		for j from i + 1 to len do
			if AreTwoWordsCompatible(wordlist[i], wordlist[j]) then
				ret1 := [ op(ret1), [i, j] ];
			fi;
		od;
		
		ret := [op(ret), ret1];
	od;
	
	ret;
end:

###################################################################################
# Verify that all triples generated from three words are legal Brinkhuis-pairs
# Assuming each triple generated from each word is permissible
# Assuming the index numbers are in ascending order
# Assuming any two words are compatible
###################################################################################

AreThreeWordsCompatible := proc(index1, index2, index3, wordlist, wordlist123, wordlist132)
	local i1, i2, i3, i, len, ret, list1, list2, list3, word1, word2, word3;
	
	word1 := wordlist[index1];
	word2 := wordlist[index2];
	word3 := wordlist[index3];
	
	list1 := [word1, word2, word3];
	list2 := [wordlist123[index1], wordlist123[index2], wordlist123[index3]];
	list3 := [wordlist132[index3], wordlist132[index3], wordlist132[index3]];
	ret := []:

	# loop through all possible triples
	# except the two generated from one word
	for i1 from 1 to 3 do
		for i2 from 1 to 3 do
			for i3 from 1 to 3 do
				# there is no need to include lists that are created by only one or two words
				if (i1 <> i2) and (i1 <> i3) and(i2 <> i3) and
						not IsTriplePermissibleForTripleWords(list1[i1], list2[i2], list3[i3], [[1,1]]) then	
					RETURN(false);
				fi;
			od;
		od;
	od;
	
	true;
end:

###################################################################################
# Given a list of words, find all triples of indexes of words 
# that are compatible with each other
###################################################################################

FindAllCompatibleTriples := proc(wordlist, compatiblepairlist)
	local ret, ret1, len, i, j, k, wordlist123, wordlist132;
	
	len := nops(wordlist);
	ret := [];
	
	wordlist123 := []:
	wordlist132 := [];
	for i from 1 to len do
		wordlist123 := [op(wordlist123), PermuteWord123(wordlist[i])];
		wordlist132 := [op(wordlist132), PermuteWord132(wordlist[i])];
	od;
		
	for i from 1 to len do
print(`FindAllCompatibleTriples: `, i, time());
		ret1 := [];
		
		for j from i + 1 to len do
print(`FindAllCompatibleTriples: `, i, j, time());

			# at least any two of the words must be compatible
			if not member([i, j], compatiblepairlist[i]) then
				next;
			fi;
			
			for k from j + 1 to len do
				# at least any two of the words must be compatible
				if not member([i, k], compatiblepairlist[i]) or
						not member([j, k], compatiblepairlist[j]) then
					next;
				fi;
			
				if AreThreeWordsCompatible(i, j, k, wordlist, wordlist123, wordlist132) then
					ret1 := [ op(ret1), [i, j, k] ];
				fi;
			od;
		od;
		
		ret := [op(ret), ret1];
	od;
	
	ret;
end:

###################################################################################
# Given a list of words, find all triples of indexes of words 
# that are compatible with each other
# Assuming reversed words are all next to each other
###################################################################################

FindAllCompatibleTriplesForWord := proc(wordindex, wordlist, compatiblepairlist, startfrom, endat)
	local ret, ret1, ret2, len, i, j, k, wordlist123, wordlist132, 
		goodwords, reversewords, normalwords, palindromes, badwords,
		isreverse1, isbadword2, isbadword3, nstart, nend, ri, rj, rk;
	
	if nargs > 5 then
		ret := args[6];
	else
#		ret := [];
	fi;
	
	goodwords := [wordindex];
	reversewords := [];
	normalwords := []:
	palindromes := []:
	badwords := []:
	
	# find all words and their reverses, and all the words whose reverses
	# are not compatible with the word
	for i from 1 to nops(compatiblepairlist[wordindex]) do
		goodwords := [op(goodwords), compatiblepairlist[wordindex][i][2]];
	od;
	len := nops(goodwords);

	for i from 1 to nops(goodwords) do
		unassign('j');
		unassign('k');

		if member(ReverseWord(wordlist[goodwords[i]]), wordlist, j) and 
				member(j, goodwords, k) then
			if i > k then
				reversewords := [op(reversewords), goodwords[i]];
			elif i < k then
				normalwords := [op(normalwords), goodwords[i]];
			else
				palindromes := [op(palindromes), goodwords[i]];
			fi;
		else
			badwords := [op(badwords), goodwords[i]];
		fi;
 	od:
	reversewords := [op(reversewords), op(palindromes)];
#print(normalwords);
#print(reversewords);
#print(palindromes);
#print(badwords);

	wordlist123 := []:
	wordlist132 := [];
	for i from 1 to nops(wordlist) do
		wordlist123 := [op(wordlist123), PermuteWord123(wordlist[i])];
		wordlist132 := [op(wordlist132), PermuteWord132(wordlist[i])];
	od;
		
print(`FindAllCompatibleTriplesForWord: total = `, len, time());
	
	if not member(startfrom, goodwords, nstart) or not	member(endat, goodwords, nend) then
		ERROR(`Illegal Start or End`);
	fi;

	for i from nstart to nend do
print(`FindAllCompatibleTriplesForWord: `, i, goodwords[i], time());
		if nops(eval(ret[goodwords[i]])) > 1 then
			next;
		fi;
	
#		ret1 := [];
		
		# we do not have to recalc everything for reversed words
		isreverse1 := member(goodwords[i], reversewords);
		
		for j from i + 1 to len do
#if i = nstart then
#print(`FindAllCompatibleTriplesForWord: `, i, goodwords[i], j, goodwords[j], time());
#fi;
			# at least any two of the words must be compatible
			if not member([goodwords[i], goodwords[j]], compatiblepairlist[goodwords[i]]) then
				next;
			elif isreverse1 then
				isbadword2 := member(j, badwords);
			fi;
			
			ret2 := [];
			for k from j + 1 to len do
				# at least any two of the words must be compatible
				if not member([goodwords[i], goodwords[k]], compatiblepairlist[goodwords[i]]) or
						not member([goodwords[j], goodwords[k]], compatiblepairlist[goodwords[j]]) then
					next;
				elif isreverse1 then
					# if each of j or k is a bad word, we have to recalc
					# otherwise, duplicate what its reverse has
					isbadword3 := member(k, badwords);
					
					if not isbadword2 and not isbadword3 then
						unassign('ri');
						unassign('rj');
						unassign('rk');
						
						member(ReverseWord(wordlist[goodwords[i]]), wordlist, ri);
						member(ReverseWord(wordlist[goodwords[j]]), wordlist, rj);
						member(ReverseWord(wordlist[goodwords[k]]), wordlist, rk);
#print(`found: `, i, j, k, ri, rj, rk, nops(ret[ri][rj]));
					
						if nops(ret[ri][rj]) > 1 and member([ri, rj, rk], ret[ri][rj]) then
#print(`in: `, i, j, k, ri, rj, rk);
							ret2 := [op(ret2), [goodwords[i], goodwords[j], goodwords[k]]];
						
							next;
						fi;
					else
						print(`been here: `, [goodwords[i], goodwords[j], goodwords[k]]);
					fi;
				fi;
			
				if AreThreeWordsCompatible(goodwords[i], goodwords[j], goodwords[k], wordlist, wordlist123, wordlist132) then
					ret2 := [op(ret2), [goodwords[i], goodwords[j], goodwords[k]]];
				fi;
			od;
			
			ret[goodwords[i]][goodwords[j]] := ret2;
#			ret1 := [op(ret1), ret2];
		od;

#		ret := [op(ret), ret1];
	od;

print(`END: `, time());	
	ret;
end:

###################################################################################
# Given a list of compatible triples, find the maximum set of words
# that are compatible with each other
###################################################################################

GetMaxTripleSet := proc(compatibletriplelist, wordlist)
	local ret, len, i, j, k, initlist, goodlist, maxlen, retword, tempword;
	
	ret := [];
	maxlen := 0;
	initlist := GetInitCompatibleWordsList(compatibletriplelist);
	len := nops(initlist);

	for i from 1 to len do
		# if the next element in the initlist is shorter than the existing result
		# exit from the loop
		if maxlen > nops(initlist[i]) then
			break;
		fi;
		
		# get the maximum set from each list
		goodlist := GetCompatibleListFromInitList(initlist[i], compatibletriplelist);

		# if the list is larger then the existing one, save it
		# if the list is of the same length as the existing one, attach it to the existing one
		if nops(goodlist[1]) > maxlen then
			maxlen := nops(goodlist[1]);
			ret := goodlist;
		elif nops(goodlist[1]) = maxlen then
			ret := [op(ret), op(goodlist)];
		fi;
	od:
	
	# change the indexes to the corrsponding words
	retword := [];
	for i from 1 to nops(ret) do
		tempword := [];
		for j from 1 to nops(ret[i]) do
			tempword := [op(tempword), wordlist[ret[i][j]]];
		od;
		
		retword := [op(retword), tempword];
	od;
	
	retword;
	ret;
end:

###################################################################################
# Given a list of compatible triples, find the lists of numbers that from the third and on,
# each number are compatible with the first two, regardless their compatibility
# with each other. It can be proven that the optimal list must be a sublist 
# of one of such lists
#
# Assume the compatible word list is in ascending order
###################################################################################

GetInitCompatibleWordsList := proc(compatibletriplelist)
	local i, j, len, ret, templist, firstnumber, secondnumber, tempword;
	
	ret := [];
	len := nops(compatibletriplelist);

	templist := [];
	firstnumber := 0;
	secondnumber := 0;
	
	for i from 1 to len do
		tempword := compatibletriplelist[i];
		
		if tempword[1] = firstnumber and tempword[2] = secondnumber then
			templist := [op(templist), tempword[3]];
		else
			if templist <> [] then
				ret := [op(ret), templist];
			fi;
			templist := tempword;
			firstnumber := tempword[1];
			secondnumber := tempword[2];
		fi;
	od;
	if templist <> [] then
		ret := [op(ret), templist];
	fi;
	
	ret := sort(ret, SortInitCompatibleList);
end:

###################################################################################
# Sort lists by their size, longest first
###################################################################################

SortInitCompatibleList := proc(list1, list2)
	if nops(list1) >= nops(list2) then
		true;
	else
		false;
	fi;
end:

###################################################################################
# Given a element from the result of  function GetInitCompatibleWordsList, which
# is also a list, find the longest compatible list from it
#
# Assume the lists are in ascending order
# the first two numbers in the list are always compatible with all the other numbers
###################################################################################

GetCompatibleListFromInitList := proc(initcompatiblelist, compatibletriples)
	local i, j, k, len, nolist, lennolist, firstnumber, secondnumber, 
		newlist, newitem, templist, maxlength, newcomplist, lennewcomplist;

	firstnumber := initcompatiblelist[1];
	secondnumber := initcompatiblelist[2];
	newlist := [];
	nolist := initcompatiblelist[3 .. nops(initcompatiblelist)];
	lennolist := nops(nolist);
	maxlength := 0;
	
	# create the first set of lists to be grown from
	# the last number is the next item to be started calculating from next time
	for i from 1 to lennolist do
		newlist := [op(newlist), [[firstnumber, secondnumber, nolist[i]], 1]];
	od;
	
	# len is the length of newlist. we keep track of it instead of 
	# using nops to improve performance
	i := 1;
	len := lennolist;
	
	while i <= len do
		newitem := newlist[i];
		templist := [];
		
		for j from newitem[2] to lennolist do
			if member(nolist[j], newitem[1]) then
				# skip the number if it is already in the list
				next;
			fi;
			
			# if the number is compatible with all the numbers in the item
			# add it in, and remember the result
			if IsNumberCompatibleWithList(nolist[j], newitem[1], compatibletriples) then
				newcomplist := sort([op(newitem[1]), nolist[j]]);
				lennewcomplist := nops(newcomplist);
				
				# make sure the resulting list has not been calculated
				# if such a list exists, ignore the last one, even if the next item index 
				# may not be the same because it will not affect the final result
				for k from len by -1 to 1 do
					if nops(newlist[k][1]) < lennewcomplist then
						k := 0;
						break;
					elif newlist[k][1] = newcomplist then
						break;
					fi;					
				od;
				
				if k < 1 then
					templist := [op(templist), [newcomplist, j + 1]];
				fi;
			fi;			
		od;
	
		if templist <> [] then
			# if the result list if not empty, we got new lists created, and they are
			# guaranteed to be at least as long as the existing ones, and longer than
			# the current one and the ones in front of it. So it is safe to remove 
			# all of them.

			newlist := [op(i + 1 .. len, newlist), op(templist)];
			len := len - i + nops(templist);
			i := 0;
			maxlength := nops(newitem[1]) + 1;
#
#			newlist := subsop(i = NULL, newlist);
#			newlist := [op(newlist), op(templist)];
#			i := i - 1;
#			len := len + nops(templist) - 1;
		else
			# if nothing new is found, then this is the (local) final result
			# but may not be the best one
			if nops(newitem[1]) < maxlength then
				newlist := subsop(i = NULL, newlist);
				len := len - 1;
				i := i - 1;
			elif nops(newitem[1]) > maxlength then
				# I do not think this part can ever be run
				newlist := newlist[i .. nops(newlist)];
				len := len - i + 1;
				i := 1;
				maxlength := nops(newitem[1]);
			fi;
		fi;
		
		i := i + 1;
	od;
	
	# remove the extra numbers (the next item indexes) at the end
	templist := [];
	for i from 1 to nops(newlist) do
		templist := [op(templist), newlist[i][1]];
	od;
	
	templist;
end:

###################################################################################
# Given a number no, find out if it is compatible with all other numbers 
# in the list nolist, that is, any triples formed using no is compatible.
# compatibletriples is the list of all compatible triples
###################################################################################

IsNumberCompatibleWithList := proc(no, nolist, compatibletriples)
	local i, j, len, newlist;

	len := nops(nolist);
	newlist := sort(nolist);
	
	for i from 1 to len do
		for j from i + 1 to len do
			if not member(sort([newlist[i], newlist[j], no]), compatibletriples) then
				RETURN(false);
			fi;
		od;
	od;
	
	true;
end:

###################################################################################
# Permute a given word with (1,2,3). Assuming 0 is not in the word
###################################################################################

PermuteWord123 := proc(word)
	local word2;
	
	word2 := word;
	word2 := ConvertLetter(word2, 1, 0):
	word2 := ConvertLetter(word2, 3, 1):
	word2 := ConvertLetter(word2, 2, 3):
	word2 := ConvertLetter(word2, 0, 2):
end:

###################################################################################
# Permute a given word with (1,3,2). Assuming 0 is not in the word
###################################################################################

PermuteWord132 := proc(word)
	local word3;

	word3 := word;
	word3 := ConvertLetter(word3, 1, 0):
	word3 := ConvertLetter(word3, 2, 1):
	word3 := ConvertLetter(word3, 3, 2):
	word3 := ConvertLetter(word3, 0, 3):
end:

###################################################################################
#
#
# OVERWRITTEN FUNCTIONS
#
#
###################################################################################


###################################################################################
# Overwritten function to determine if a word has a square in it
#
# word is a list of letters
# patterns is assumed to be [[1,1]], thus ignored
# solidbadwords is ignored
###################################################################################

IsBadWord := proc(word, patterns, solidbadwords)
	local len, i, ii, j, replacelist, tempword;

	len := nops(word);
	
	for i from 1 to len do
		for j from i to (len + i - 1) / 2 do
			if word[i .. j] = word[j + 1 .. 2 * j - i + 1] then
				RETURN(true);
			fi;
		od;
	od;
	
	false:
end:

###################################################################################
# A variation of IsTriplePermissible
#
# The function is to be called in functions where the set of word2 and word3 will 
# be images of the set of word1 under certain permutations. 
# Hence, when we try to verify the result,
# we can make word1 always be the first part in the combined string
# because all other cases will be encountered as a image of a permutation sooner or later,
# which does not affect the compatibility
###################################################################################

IsTriplePermissibleForCompatibleWords := proc(word1, word2, word3, patterns)
	local i, j, goodlist, lenlist, conversion, result;

	goodlist := [[1, 2, 3], [1, 3, 2], [1, 2, 1], [1, 3, 1]]:
	lenlist := nops(goodlist);
	conversion := [[[1], [word1]], [[2], [word2]], [[3], [word3]]];
	
	for i from 1 to lenlist do
		result := ConvertString(goodlist[i], conversion)[1];

		if IsBadWord(result, patterns, []) then
			RETURN(false);
		fi;
	od;
	
	TestExtraCases(word1, word2, word3, patterns);
end:

###################################################################################
# A variation of IsBadWord
#
###################################################################################

IsBadWordForTripleWords := proc(word, sectionlength)
	local len, i, j, k;

	len := 3 * sectionlength;
	
	# a word has three sections of equal lengths, of which any two are compatible
	# to find out the whole word is bad or not, we only need to check 
	# if two consecutive subwords are equal where
	# the start of the first is in the first section
	# and the end of the second is in the third section
	# i.e. i <= sectionlength, and 2 * j - i + 1 >= 2 * sectionlength + 1
	for i from 1 to sectionlength do
		for j from sectionlength + iquo(i + 1, 2) to (len + i - 1) / 2 do
			if word[i .. j] = word[j + 1 ..  2 * j - i + 1] then
				RETURN(true);
			fi;
		od;
	od;
	
	false:
end:

###################################################################################
# A variation of IsTriplePermissible
#
# The function is to be called in functions where the set of word2 and word3 will 
# be images of the set of word1 under certain permutations. 
# Hence, when we try to verify the result,
# we can make word1 always be the first part in the combined string
# because all other cases will be encountered as a image of a permutation sooner or later,
# which does not affect the compatibility
###################################################################################

IsTriplePermissibleForTripleWords := proc(word1, word2, word3, patterns)
	local i, j, goodlist, lenlist, conversion, result, lenword;

	lenword := nops(word1);
		
	if IsBadWordForTripleWords([op(word1), op(word2), op(word3)], lenword) then
		RETURN(false);
	elif IsBadWordForTripleWords([op(word1), op(word3), op(word2)], lenword) then
		RETURN(false);
	elif IsBadWordForTripleWords([op(word1), op(word2), op(word1)], lenword) then
		RETURN(false);
	elif IsBadWordForTripleWords([op(word1), op(word3), op(word1)], lenword) then
		RETURN(false);
	fi;

	TestExtraCases(word1, word2, word3, patterns);
end:


###################################################################################
#
#
# HEAD AND TAILS
#
#
###################################################################################

###################################################################################
# Given a head and a tail, find a new tail that extends tail by one letter at the
# beginning so that both [new_tail, PermuteWord123(head)] 
# and [new_tail, PermuteWord132(head)] are good words
###################################################################################

FindNewTail := proc(head, tail)
	local i, j, retgood, retbad, no, th, head1, head2, len, newtail, result;
	
	retgood := [];
	retbad := [];
	
	head1 := PermuteWord123(head);
	head2 := PermuteWord132(head);
	
	if nops(tail) = 0 then
		newtail := {1, 2, 3} minus {head1[1], head2[1]};

		for i from 1 to nops(newtail) do
			retgood := [op(retgood), [head, [newtail[i]]]];
		od;
	else
		th := tail[1];
		
		# we only need to test the two numbers 
		# that are not the same as the beginning of tail
		no := (th mod 3) + 1;
		result := TestLetterForTail(no, tail, head1, head2);
		if result = 1 then
			retgood := [op(retgood), [head, [no, op(tail)]]];
		elif result = -1 then
			retbad := [op(retbad), [head, [no, op(tail)]]];
		fi;
		
		no := ((th + 1) mod 3) + 1;
		result := TestLetterForTail(no, tail, head1, head2);
		if result = 1 then
			retgood := [op(retgood), [head, [no, op(tail)]]];
		elif result = -1 then
			retbad := [op(retbad), [head, [no, op(tail)]]];
		fi;
	fi;
	
	retgood, retbad;
end:

###################################################################################
# test if a number no is good if added to the beginning of a tail
#
# return 1 if good, 
# return 0 if the newtail itself is bad, 
# return -1 if the newtail is good, but no compatible with the heads
###################################################################################

TestLetterForTail := proc(no, tail, head1, head2)
	local i, len, newtail;

	newtail := [no, op(tail)];
		
	if IsBadWord(newtail) then
#PrintResult(newtail, ` is not allowed because itself is bad`);
		RETURN(0):
	else
		if IsBadWord([op(newtail), op(head1)]) or
				IsBadWord([op(newtail), op(head2)]) then
#PrintResult(newtail, ` is not allowed because it is not compatible with heads`);
			RETURN(-1);
		fi;

		RETURN(1);			
	fi;
end:

###################################################################################
# Given a head and a tail, find a new head that extends head by one letter at the
# end so that both [PermuteWord123(tail), new_head] 
# and [PermuteWord132(tail), new_head] are good words
###################################################################################

FindNewHead := proc(head, tail)
	local i, j, retgood, retbad, no, th, tail1, tail2, len, newhead, result;
	
	retgood := [];
	retbad := [];
	
	tail1 := [op(tail1), PermuteWord123(taillist[i])];
	tail2 := [op(tail2), PermuteWord132(taillist[i])];
	
	if nops(head) = 0 then
		newhead := {1, 2, 3} minus {tail1[1], tail2[1]};

		for i from 1 to nops(newhead) do
			retgood := [op(retgood), [head, [newhead[i]]]];
		od;
	else
		th := head[nops(head)];
		
		# we only need to test the two numbers 
		# that are not the same as the end of head
		no := (th mod 3) + 1;
		result := TestLetterForHead(no, head, tail1, tail2);
		if result = 1 then
			retgood := [op(retgood), [[op(head), no], tail]];
		elif result = -1 then
			retbad := [op(retbad), [[op(head), no], tail]];
		fi;
		
		no := ((th + 1) mod 3) + 1;
		result := TestLetterForHead(no, head, tail1, tail2);
		if result = 1 then
			retgood := [op(retgood), [[op(head), no], tail]];
		elif result = -1 then
			retbad := [op(retbad), [[op(head), no], tail]];
		fi;
	fi;
	
	retgood, retbad;
end:

###################################################################################
# test if a number no is good if added to the end of a head
#
# return 1 if good, 
# return 0 if the newhead itself is bad, 
# return -1 if the newhead is good, but no compatible with the tails
###################################################################################

TestLetterForHead := proc(no, head, tail1, tail2)
	local i, len, newhead;

	newhead := [op(head), no];
		
	if IsBadWord(newhead) then
#PrintResult(newhead, ` is not allowed because itself is bad`);
		RETURN(0):
	else
		if IsBadWord([op(tail1), op(newhead)]) or
				IsBadWord([op(tail2), op(newhead)]) then
#PrintResult(newhead, ` is not allowed because it is not compatible with heads`);
			RETURN(-1);
		fi;

		RETURN(1);			
	fi;
end:

###################################################################################
# Find all the comibinations of good and bad heads and tails up to length letters.
###################################################################################

FindHeadsAndTails := proc(length)
	local i, j, k, retgood, retbad, tempgood, tempbad, temphead, temptail;
	
	retgood := [[[1,2], [2,1]]];
	retbad := [];
	
	for i from 3 to length do
		tempgood := retgood;
		retgood := [];
		
		for j from 1 to nops(tempgood) do
			temphead := FindNewHead(op(tempgood[j]));
			retbad := [op(retbad), op(temphead[2])];
			
			for k from 1 to nops(temphead[1]) do
				temptail := FindNewTail(op(temphead[1][k]));
				retgood := [op(retgood), op(temptail[1])];
				retbad := [op(retbad), op(temptail[2])];
			od;
		od;
	od;
	
	retgood, retbad;
end:





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

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

GetMaxTripleSet := proc(compatibletriplelist, wordlist, startfrom, endto)
	local ret, len, i, j, k, initlist, compatiblelist, maxlen, retword, tempword;
	
	ret := [];
	maxlen := 0;
	initlist := GetInitCompatibleWordsList(compatibletriplelist, startfrom, endto);
	len := nops(initlist);

	for i from 1 to len do
		# if the next element in the initlist is shorter than the existing result
		# exit from the loop
		if maxlen > initlist[i][3] + 2 then
			break;
		fi;
		
		compatiblelist := [];
		for j from 1 to nops(compatibletriplelist[initlist[i][1]][initlist[i][2]]) do
			compatiblelist := [op(compatiblelist), compatibletriplelist[initlist[i][1]][initlist[i][2]][j][3]];
		od;
		
		# get the maximum set from each list
		ret := FindMaxCompatibleList([initlist[i][1], initlist[i][2]], compatiblelist, compatibletriplelist, ret);
		maxlen := nops(ret[1]);
print(`GetMaxTripleSet: `, i, len, maxlen, time());		
	od:
	
	# change the indexes to the corrsponding words
	retword := [];
	for i from 1 to nops(ret) do
		tempword := [];
		for j from 1 to nops(ret[i]) do
			tempword := [op(tempword), wordlist[ret[i][j]]];
		od;
		
		retword := [op(retword), tempword];
	od;
	
	retword;
	
	ret;
end:

###################################################################################
# Given a list of compatible triples, find the lists of numbers that from the third and on,
# each number are compatible with the first two, regardless their compatibility
# with each other. It can be proven that the optimal list must be a sublist 
# of one of such lists
#
# Assume the compatible word list is in ascending order
###################################################################################

GetInitCompatibleWordsList := proc(compatibletriplelist, startfrom, endto)
	local i, j, len, ret, templist, firstnumber, secondnumber, tempword;
	
	ret := [];

	templist := [];
	firstnumber := 0;
	secondnumber := 0;
	
	for i from startfrom to endto do
		for j from i + 1 to endto do
			if type(compatibletriplelist[i][j], list) then
				ret := [op(ret), [i, j, nops(compatibletriplelist[i][j])]];
			fi;
		od;
	od;
	
	ret := sort(ret, SortInitCompatibleList);
end:

###################################################################################
# Sort lists by the third elemenent descending, then the first and second, ascending
###################################################################################

SortInitCompatibleList := proc(list1, list2)
	if list1[3] > list2[3] then
		true;
	elif list1[3] < list2[3] then
		false;
	elif list1[1] < list2[1] then
		true;
	elif list1[1] > list2[1] then
		false;
	elif list1[2] < list2[2] then
		true;
	else
		false;
	fi;
end:

FindMaxCompatibleList := proc(goodlist, compatiblelist, compatibletriplelist, maxlist)
	local i, j, len1, len2, tempgoodlist, tempcomplist, 
		maxlength, tempset, tempall, lenmax, ret;
	
	len1 := nops(goodlist);
	len2 := nops(compatiblelist);
	lenmax := nops(maxlist);
	ret := maxlist;
	
	if lenmax > 0 then
		maxlength := nops(maxlist[1]);
	else
		maxlength := 0;
	fi;
	
#	if len1 + len2 < maxlength then
#		RETURN(maxlist);
#	fi;
		
	if len2 = 0 then
		if len1 + len2 > maxlength then
PrintResult(`FindMaxCompatibleList: found new max result: `, [op(goodlist), op(compatiblelist)], time());
			RETURN([[op(goodlist), op(compatiblelist)]]);
		elif len1 + len2 = maxlength then
			for i from 1 to lenmax do
				ret := {op(goodlist), op(compatiblelist)};
				
				if ret = {op(maxlist[i])} then
					break;
				fi;
			od;
			
			if i <= lenmax then
				RETURN(maxlist);
			else
#PrintResult(`FindMaxCompatibleList: found another result: `, [op(goodlist), op(compatiblelist)], time());
				RETURN([op(maxlist), [op(goodlist), op(compatiblelist)]]);
			fi;
		else
			RETURN(maxlist);
		fi;
	fi;
	
	for i from 1 to len2 do
		tempset := [];
		for j from 1 to nops(ret) do
			tempset := [op(tempset), {op(ret[j])}];
		od;
		
		tempgoodlist := [op(goodlist), compatiblelist[i]];
		tempcomplist := GetNextCompatibleList(tempgoodlist, 
			subsop(i = NULL, compatiblelist), compatibletriplelist);
		
		tempall := {op(tempgoodlist), op(tempcomplist)};
		if nops(tempall) >= maxlength then
			for j from 1 to lenmax do
				if tempall minus tempset[j] = {} then
					break;
				fi;
			od;
			
			if j <= lenmax then
				next;
			fi;
			
			ret := FindMaxCompatibleList(tempgoodlist, tempcomplist, compatibletriplelist, ret);
			maxlength := nops(ret[1]);
		fi;
	od;
	
	ret;
end:

GetNextCompatibleList := proc(goodlist, compatiblelist, compatibletriplelist)
	local i, len, ret, ret1;

	len := nops(compatiblelist);
	ret := [];
	
	for i from 1 to len do
		if IsNumberCompatibleWithList(compatiblelist[i], goodlist, compatibletriplelist) then
			ret := [op(ret), compatiblelist[i]];
		fi;
	od;

	ret;
end:

###################################################################################
# Given a number no, find out if it is compatible with all other numbers 
# in the list nolist, that is, any triples formed using no is compatible.
# compatibletriplelist is the list of all compatible triples
###################################################################################

IsNumberCompatibleWithList := proc(no, nolist, compatibletriplelist)
	local i, j, len, newlist, temp;

	len := nops(nolist);
	newlist := sort(nolist);
	
	for i from 1 to len do
		for j from i + 1 to len do
			temp := sort([newlist[i], newlist[j], no]);
			
			if not type(compatibletriplelist[temp[1]][temp[2]], list) then
				RETURN(false);
			elif not member(temp, compatibletriplelist[temp[1]][temp[2]]) then
				RETURN(false);
			fi;
		od;
	od;
	
	true;
end:

AddValueToEachElement := proc(datalist, value)
	local i, j, ret, tt;
	
	ret := [];
	
	for i from 1 to nops(datalist) do
		if type(datalist[i], list) or type(datalist[i], set) then
			tt := AddValueToEachElement(datalist[i], value);
		else
			tt := datalist[i] + value;
		fi;
		
		ret := [op(ret), tt];
	od;
	
	ret;
end: