##############################################################################
# we define a Chomp position as a ordered list, 
# e.g. a three-row position [a1, a2, a3] looks like
# xxx					the first row  = a1
# xxxxxx				the second row = a1 + a2
# xxxxxxxxxxx			the third row = a1 + a2 + a3
#
# A nim value is the value of the Sprague-Grundy function
#
# InitNim is a list of nim values that InitNim[k][h] is a list
# of nim values with k rows and the value of the first row equals to h
#
# this program returns the nim values (Grundy function) for all positions
##############################################################################

print(`Type \'Help()\' for help`);

InitNim := [[[[1+a[1], a[1]]]]]:

#MAX_LAST_ROW_SIZE := 250:
MAX_POLY_SIZE := 50000000:
MIN_OCCURENCE := 3:
Symb_Var := `x`: #'Symb_Var':
Symb_Nim := `y`:
CHECK_VALUE := 1000000:
PRINT_RESULT := true:

NON_FIXED_ROWS := 1:

readlib(mtaylor):


##############################################################################
# The help function
# Usage: Help for the help menu
#        Help(functionname) for the specific help for the function
##############################################################################

Help := proc()
	if args = NULL then
		printf(`%s\n`, `This program calculates nim values (Sprague-Grundy function) for all Chomp positions`);
		printf(`%s\n`, `The main function is Main`);
		printf(`%s\n`, ``);
		printf(`%s\n`, `The nim value for a Chomp position is an ordered list`);
		printf(`%s\n`, `the number of operand is one more than the number of rows of the position`);
		printf(`%s\n`, `the first operand is the number of pieces in the first row`);
		printf(`%s\n`, `followed by the numbers of excessive pieces in `);
		printf(`%s\n`, `the corresponding rows, while the last number is the nim value for the position. `);
		printf(`%s\n`, `For example [2, 0, 4, 8] represents a position `);
		printf(`%s\n`, `with 3 rows (2, 2, 6 pieces) and 8 as its nim value`);
		printf(`%s\n`, ``);
		printf(`%s\n`, `Help is available for the following functions:`);
		printf(`%s\n`, `Constants, Main, FindPositionNimValue, FindNextMoveWithNim, PlayChomp`);
		printf(`%s\n`, ``);
		printf(`%s\n`, `Usage: Help(funtionname) for the help of the function functionname`);
	elif nops(args) = 1 then
		if op(1, [args]) = Main then
			printf(`%s\n`, `Main function takes four perameters: T, max_rows, max_cols, and a`);
			printf(`%s\n`, `T is the precalculated P-positions with predefined format.`);
			printf(`%s\n`, `        Please use InitNim or results from previous calculations`);
			printf(`%s\n`, `max_rows is the maximum number of rows that the Chomp positions are to have`);
			printf(`%s\n`, `max_cols is the maximum number of columns that the Chomp positions the top `);
			printf(`%s\n`, `        (max_rows - 2) rows are to have`);
			printf(`%s\n`, `a is a symbolic integer that the formulas will use.`);
			printf(`%s\n`, ``);
			printf(`%s\n`, `In other words, we are currently calculating P-positions with `);
			printf(`%s\n`, `the following pattern:`);
			printf(`%s\n`, ``);
			printf(`%s\n`, `*                             `);
			printf(`%s\n`, `*                             `);
			printf(`%s\n`, `**                            `);
			printf(`%s\n`, `******                        `);
			printf(`%s\n`, `******                        `);
			printf(`%s\n`, `*********************         `);
			printf(`%s\n`, `******************************`);
			printf(`%s\n`, ``);
			printf(`%s\n`, `The format for the input T and output will be the same:`);
			printf(`%s\n`, `It is an ordered list whose operands are ordered lists in which`);
			printf(`%s\n`, `the i-th operand is the list of P-positions with the i rows,`);
			printf(`%s\n`, `and the j-the operand of this list is the list of `);
			printf(`%s\n`, `all P-positions with i rows and the value of the first row is j.`);
			printf(`%s\n`, ``);
			printf(`%s\n`, `For example, there is a pre-defined list InitNim = [[[[1]]], [[[1 + a, 1]]]]`);
			printf(`%s\n`, `which means, all the P-positions with one row is [1],`);
			printf(`%s\n`, `all the P-positions with two rows has the form of [1 + a, 1], `);
			printf(`%s\n`, `where a is a symbolic integer greater than or equal to 0`);
			printf(`%s\n`, ``);
			printf(`%s\n`, `Users can plug in a pre-calculated list T as the starting point.`);
		elif op(1, [args]) = FindPositionNimValue then
			printf(`%s\n`, `The function takes two perameters: POS and a`);
			printf(`%s\n`, `where POS is the Chomp position whose nim value is to be found`);
			printf(`%s\n`, `a is a symbolic integer that the result set (i.e. InitNim) is using`);
		elif op(1, [args]) = FindNextMoveWithNim then
			printf(`%s\n`, `The function takes three perameters: POS, val and a`);
			printf(`%s\n`, `where POS is the Chomp position whose next move with nim value val is to be found`);
			printf(`%s\n`, `a is a symbolic integer that the result set (i.e. InitNim) is using`);
		elif op(1, [args]) = PlayChomp then
			printf(`%s\n`, `The function takes one perameter: POS`);
			printf(`%s\n`, `where POS is the Chomp position to start for the game`);
			printf(`%s\n`, `Human always has the first move`);
			printf(`%s\n`, `Following the instruction to play the game. Whoever claims 1, 1 loses,`);
			printf(`%s\n`, `or in other words, whoever leaves the game with only one piece wins.`);
		elif op(1, [args]) = Constants then
			printf(`%s\n`, `The following are predefined constants used in the program:`);
			printf(`%s\n`, ``);
			printf(`%s%a;\n`, `InitNim = `, InitNim);
			printf(`%s\n`, ``);
			printf(`%s\n`, `InitNim is the default nim value list. See help for the Main function for detail`);
			printf(`%s\n`, ``);
			printf(`%s%d;\n`, `MIN_OCCURENCE = `, MIN_OCCURENCE);
			printf(`%s\n`, ``);
			printf(`%s\n`, `MIN_OCCURENCE is the minimal occurance of a pattern to be `);
			printf(`%s\n`, `considered as a formula. The default value of MIN_OCCURENCE `);
			printf(`%s\n`, `is 4, which means certain type of value must happen at least three `);
			printf(`%s\n`, `times to be considered as a possible formula,  e.g. 24, 26, 28, 30 is `);
			printf(`%s\n`, `enough to be considered as 24 + 2a when MIN_OCCURENCE = 4`);
			printf(`%s%a;\n`, `PRINT_RESULT = `, PRINT_RESULT);
			printf(`%s\n`, ``);
			printf(`%s\n`, `PRINT_RESULT is used to print redundant messages.`);
			printf(`%s\n`, `Redundant messages is printed if PRINT_RESULT is true and printlevel > 1.`);
			printf(`%s\n`, `The default value is false.`);
			printf(`%s\n`, ``);
		else
			printf(`%s\n`, `Unknown function name. Available functions are: `);
			printf(`%s\n`, `        Constants, Main, FindSymbolicWinningMove, FindWinningMove`);
		fi:
	else
		printf(`%s\n`, `Usage: Help(funtionname) for the help of the function functionname`);
	fi:
end:

##############################################################################
# Display all pre-defined constants
# Reuse the function body in Help(Constants)
##############################################################################

Constants := proc()
	Help(Constants);
end:

##############################################################################
# Main Function
#
# InitNimos are the known nim values for Chomp positions
# max_rows is the max number of rows to calculate
# max_cols is the max number of cols for each row
# a is the symbolic variable
#
# returns Chomp positions with their nim values
##############################################################################

Main := proc(InitNimos, max_rows, max_cols, max_top, a)
	local T, i, j, toplist, top, len;
	
	T := InitNimos;
	
	for i from 2 to max_rows do
		# Get all possible top rows	
		toplist := sort(convert(ExpandBlockByMaxValue([], max_cols, 
				[max_cols], i - NON_FIXED_ROWS, a), list), SortChomp);
		len := nops(toplist);		

		# Find the nim values
		for j from 1 to len do
			if toplist[j][1] <> 0 and SortChomp([op(1..i - NON_FIXED_ROWS, toplist[j])], max_top) then
				print(cat(`We are calculating positions with top rows: `, 
					[op(1..i - NON_FIXED_ROWS, toplist[j])]));

				T := GetNimForNextFixedTop(T, [op(1..i - NON_FIXED_ROWS, toplist[j])], 
						max_cols, a);
			fi;
		od:
	od:
	
	T;
end:

##############################################################################
# Get the nim value for the next Chomp position whose top rows are FixedTop
# calculate up to max_cols number of columns
#
# T is the list of known nim values for positions
# FixedTop is the top rows for the positions to be calculated
# max_cols is the max number of columns to be calculated
# a is the symbolic variable
#
# returns all nim values for the positions with the FixedTop
##############################################################################

GetNimForNextFixedTop := proc(T, FixedTop, max_cols, a)
	local i, j, k, varlist, numerf, denomf, maxdeg, degx, degy, list1, itemf, coefff,
		NimTable, POS, f, f1, listf, fsize, genf, jj, POS1, len, bIntersect, lastremovedelement, 
		convf, AllNList, NList, rows, NimTable1, AllSolidPositions, NumOfSymbPos, status6;

	AllSolidPositions := [];
	if nops(T) >= nops(FixedTop) + NON_FIXED_ROWS then
		if nops(T[nops(FixedTop) + NON_FIXED_ROWS]) >= FixedTop[1] then
			NimTable := T[nops(FixedTop) + NON_FIXED_ROWS][FixedTop[1]];	
			len := nops(NimTable);
			
			for i from 1 to len do
				POS := NimTable[i];
				
				if [op(1..nops(FixedTop), POS)] = FixedTop then
					AllSolidPositions := [op(AllSolidPositions), POS];
				fi:
			od:
		fi:
	fi:
	
	NimTable := T;
	genf := BuildGeneratingFunctionForFixedTop(NimTable, FixedTop, max_cols, a);
	AllNList := genf[2];
	genf := genf[1];
	varlist := {};
	rows := nops(FixedTop) + NON_FIXED_ROWS;
	NumOfSymbPos := 0;
	
	for i from 1 to NON_FIXED_ROWS do
		varlist := varlist union {Symb_Var[nops(FixedTop) + i]};
	od;
	varlist := varlist union {Symb_Nim};

	# coefff is the coeff for the next Chomp position with positive coefficient
	coefff := 1;
	while coefff >=1 do
		numerf := numer(genf);
		denomf := denom(genf);
	
		# we only need to calculate up to maxdeg degrees
		maxdeg := degree(numerf) + degree(denomf) + 2;

		listf := GetItemsFromPolynomial(genf);
		fsize := 1:
#if FixedTop = [4, 0] and nops(AllSolidPositions) > 60 then
#print(` before mtaylor: number of items = `, nops(genf), nops(listf));
#fi:

		if nops(listf) > MAX_POLY_SIZE then
			f := 0;
			
			while nops(listf) > 0 do
#				f1 := sum(listf[ii], ii=1..min(MAX_POLY_SIZE, nops(listf)));

				f1 := 0;
				
				for i from 1 to min(MAX_POLY_SIZE, nops(listf)) do
					f1 := f1 + listf[i];
				od:

#if FixedTop = [4, 0] then
#print(` mtayloring `, f1, maxdeg, nops(f1));
#fi:

				f[fsize] := 0;
				status6 := 1000000;
				while status6 > 10000 do
					gc();
					
					status6 := status[6];
				od:				

				f[fsize] := mtaylor(f1, varlist, maxdeg);
				fsize := fsize + 1;
				
#if FixedTop = [4, 0] then
#print(` mtaylored `, f1);
#fi:

#				f := f + f1;
				
				listf := [op(MAX_POLY_SIZE + 1..nops(listf), listf)];
			od:

			fsize := fsize - 1;
		else
			f[fsize] := 0;
#			status6 := 1000000;
#			while status6 > 10000 do
#				gc();
#				
#				status6 := status[6];
#			od:				

			f[fsize] := mtaylor(genf, varlist, maxdeg);
		fi:

		#### !!!! CORRECT ONLY WHEN NON_FIXED_ROWS = 1 !!!!
		if NON_FIXED_ROWS <> 1 then
			ERROR(`This program cannot handle more than one non-fixed rows`);
		fi;
		
#		degx := min(degree(f, Symb_Var[nops(FixedTop) + 1]), max_cols - sum(FixedTop[ii], ii=1..nops(FixedTop)));
		degx := max(seq(degree(f[i], Symb_Var[nops(FixedTop) + 1]), i=1..fsize));
		degy := max(seq(degree(f[i], Symb_Nim), i=1..fsize));
		coefff := 0;
		
		# Find the next position and its nim value
		for i from 0 to degx do
#			itemf := sum(coeff(f[ii], Symb_Var[nops(FixedTop) + 1], i), ii=1..fsize);

			itemf := 0;
			for j from 1 to fsize do
				itemf := itemf + coeff(f[j], Symb_Var[nops(FixedTop) + 1], i);
			od;

			for j from 0 to degy do
				coefff := coeff(itemf, Symb_Nim, j):

				if coefff >= 1 then
					break;
				fi:
			od:
				
			if coefff >= 1 then
				break;
			fi:
		od:

		for k from 1 to fsize do
			f[k] := 0;
		od;

		if coefff >= 1 then
			# If we find one, we add it to our list 
			# and check the validity of the previous formulae
			
			POS := [op(FixedTop), i, j];
			AllSolidPositions := [op(AllSolidPositions), POS]; # remember all solid positions
			NimTable := AddNewPPositionIntoResult(NimTable, [POS], 
				nops(FixedTop)+1, FixedTop[1], AllNList, a);

#if FixedTop = [4, 0] and POS[rows]>=152 then
#RETURN(NimTable[2]);
#fi:

			POS := NimTable[1][1];
			NimTable := NimTable[2];
			convf := ConvertPositionToGeneratingFunctionWithFixedTop([op(1..nops(POS)-1, POS)], POS, a):
			NList := convf[2]; 
			
			len := nops(NimTable[nops(POS)-1][POS[1]]);
			bIntersect := false;
			lastremovedelement := len;
			
			if degree(POS[rows - NON_FIXED_ROWS + 1]) > 0 then
				NumOfSymbPos := NumOfSymbPos + 1;
			fi:
			
			for j from 1 to len do
				POS1 := NimTable[nops(POS)-1][POS[1]][j];
				
				if POS <> POS1 
					and (ElementsDisagree(POS, [POS1], a)
					or ElementsIntersect(POS1, NList, a)) then

					# remove the element because the pattern we guessed conflict with 
					# the newly found position
					if type(NimTable[rows][FixedTop[1]][j][rows - NON_FIXED_ROWS + 1], integer) then
						# Solid position cannot be removed because if the new nim value
						# is less than the old value, then what made the old value incorrect?
						# Probably it was because of an incorrect pattern, but it shall have\
						# be corrected before the current position is calculated.
						# On the other hand, if the new value is bigger, then the 
						# NList must have missed something the make the old value
						# so small.
						# So if we are removing a solid position instead of a pattern,
						# something must be wrong
						ERROR(cat(`Removing a solid position`, POS1, ` by `, POS, `. Something is wrong.`));
					fi:

					lastremovedelement := min(lastremovedelement, j);
					
					if NumOfSymbPos > 5 then
						# Since we have found a significant number of patterns, we only need to
						# remove the pattern and add in MIN_OCCURENCE items in the pattern
						# to speed up the program (by a little) 

PrintResult(cat(`!!!!! GetNimForNextFixedTop: element `, POS1, ` REMOVED by `, POS)); #, NimTable)):
	
						NimTable := [ op( 1..rows - 1, NimTable ),
							[ op( 1..FixedTop[1] - 1, NimTable[rows] ),
							sort( [op( 1..j - 1, NimTable[rows][FixedTop[1]] ),
							seq( subs(a[rows - NON_FIXED_ROWS + 1] = jj, NimTable[rows][FixedTop[1]][j]), jj = 0..MIN_OCCURENCE - 1 ),
							op( j + 1..nops(NimTable[rows][FixedTop[1]]), NimTable[rows][FixedTop[1]] ) ], SortChomp),
							op( FixedTop[1] + 1..nops(NimTable[rows]), NimTable[rows] ) ],
							op( rows + 1..nops(NimTable), NimTable ) ];
						
					else
PrintResult(cat(`!!!!! GetNimForNextFixedTop All Symbolic Positions are removed: element `, POS1, ` REMOVED by `, POS)); #, NimTable)):

						# Since all the patterns are compatible with the one being removed
						# we might as well simply remove all the patterns with the FixedTop

						NimTable1 := NimTable[rows][FixedTop[1]];

						len := nops(NimTable1);
						i := 1;
						while i <= len do
							POS1 := NimTable1[i];

							if [op(1..rows - NON_FIXED_ROWS, POS1)] = FixedTop then
								NimTable1 := [op( 1..i - 1, NimTable1 ),
									op( i + 1..len, NimTable1 ) ];

								len := len - 1:
	#PrintResult(cat(`Relevant symbolic result is removed: `, POS1, NimTable1));
							else
								i := i + 1;		
							fi:
						od:
							
						# restore all solid positions, at least they are all correct
						NimTable1 := [op(NimTable1), op(AllSolidPositions)]:
										
						NimTable := [ op( 1..rows - 1, NimTable ),
							[ op( 1..FixedTop[1] - 1, NimTable[rows] ),
							sort( NimTable1, SortChomp ),
							op( FixedTop[1] + 1..nops(NimTable[rows]), NimTable[rows] ) ],
							op( rows + 1..nops(NimTable), NimTable ) ];
					fi:

					bIntersect := true;
				fi:
			od:
			
			if bIntersect then
				# Since we just removed a pattern, we need to recalculate 
				# the generation function
				genf := BuildGeneratingFunctionForFixedTop(NimTable, FixedTop, max_cols, a);
				AllNList := genf[2];
				genf := genf[1];

				list1 := NimTable[nops(FixedTop) + NON_FIXED_ROWS][FixedTop[1]];
				len := nops(list1);
				i := 1:	#lastremovedelement;
				
				while i <= len do
					POS1 := list1[i];
					
					if degree(POS1[nops(FixedTop)+1]) = 0 then
						POS1 := FindPattern(POS1, list1, AllNList, false, a);
						
						if degree(POS1[1][nops(FixedTop)+1]) > 0 then
							NimTable := [seq(NimTable[k], k = 1..rows - 1), 
								[seq(NimTable[rows][k], k = 1..FixedTop[1]-1), POS1[2], 
								seq(NimTable[rows][k], k = FixedTop[1]+1..nops(NimTable[rows]))],
								seq(NimTable[k], k = rows + 1..nops(NimTable) )];

							convf := ConvertPositionToGeneratingFunctionWithFixedTop(
								[op(1..nops(POS1[1])-1, POS1[1])], POS1[1], a):
							NList := convf[2]; 

							AllNList := [op(AllNList), op(NList)];
							genf := genf - convf[1];

							list1 := NimTable[nops(FixedTop)+1][FixedTop[1]];
							len := nops(list1);
						fi:
					fi:
					
					i := i + 1;
				od:
			else
				# Add the newly found nim value into the generating funciton
				AllNList := [op(AllNList), op(NList)];
				genf := genf - convf[1];
			fi:
		fi:
	od:	
	
	NimTable;
end:

##############################################################################
# Add new positions into the nim value list
# also try to check if there are patterns forming
#
# returns a list of two elements
# the first one is the list of newly added nim values
# the second is the list of the new nim value list
#
# NimTable is the list of nim values
# value is a list of nim values to be added
# rows is the number of rows the new positions have
# cols is the value of the first row
# NList is the list of "forbidden nim values"s
# a is a list of symbolic variables
##############################################################################

AddNewPPositionIntoResult := proc(NimTable, value, rows, cols, NList, a)
	local i, j, k, POSAdd, list1, NimTable1, returnlist;
#print(`AddNewPPositionIntoResult ENTERING `, value):
	
	if nops(NimTable) < rows then
PrintResult(cat(`AddNewPPositionIntoResult new FIRST ROW item added `, op(value))):
		returnlist := sort(value, SortChomp);
		NimTable1 := [ returnlist ];
	elif nops(NimTable) >= rows and nops(NimTable[rows]) < cols then
PrintResult(cat(`AddNewPPositionIntoResult new FIRST COL item added `, op(value))):
		returnlist := sort(value, SortChomp):
		NimTable1 := [ seq(NimTable[rows][k], k=1..cols - 1), returnlist ];
	else
		returnlist := [];
		list1 := NimTable[rows][cols];
		
		for j from 1 to nops(value) do
			POSAdd := value[j];	# the tested one

			list1 := FindPattern(POSAdd, list1, NList, true, a);
			returnlist := [op(returnlist), list1[1]];
			list1 := list1[2];
		od:
#print(`AddNewPPositionIntoResult list1 = `, list1); 
		
		NimTable1 := [seq(NimTable[rows][k], k = 1..cols -1), list1, 
			seq(NimTable[rows][k], k = cols + 1..nops(NimTable[rows]))];
	fi:			

#print(`AddNewPPositionIntoResult returnlist = `, returnlist, `result = `, [seq(NimTable[k], k = 1..rows - 1), 
#		NimTable1, seq(NimTable[k], k = rows + 1..nops(NimTable) )]):
	[returnlist, [seq(NimTable[k], k = 1..rows - 1), 
		NimTable1, seq(NimTable[k], k = rows + 1..nops(NimTable) )]];
end:

##############################################################################
# Given a solid nim value, try to find out if there is a pattern forming in 
# the nim value PList. A pattern is when certain criteria happens at least
# MIN_OCCURENCE times, and the pattern does not conflict with any
# of the "forbidden nim values"s in NList.
#
# returns a list of two elements. The first element is the newly added
# Chomp position (may be symbolic if a pattern is found), the second
# is the new PList
#
# POSAdd is the newly found position and its nim value
# PList is the list of all nim values
# NList is the list of all "forbidden nim values"
# a is a list of symbolic variables
##############################################################################

FindPattern := proc(POSAdd, PList, NList, bAddPOSIfNoPatternFound, a)
	local len, gap, firstoccurenceat, occurence, firstindex, i, j, POS1,
		POS, rows, lastgap, patternlist, resultlist, sol, redudantlist;
	
	len := nops(PList);
	rows := nops(POSAdd);
	
	lastgap := 0;
	patternlist := [];
	firstoccurenceat := len;

	# Test all the possibilities until a pattern is found
	while true do	# get all the possible patterns available
		firstindex := POSAdd[rows - 1];
		gap := 0:
		occurence := 1:

		# Try to find a pattern in the positions, regardless of their nim values
		for i from firstoccurenceat by -1 to 1 do
			POS := PList[i]:
			
			if not type(POS[rows - 1], integer) then
				next:
			fi:
			
			if gap = 0 then
				# find the first possible position that might generate a symbolic sequence
				# firstindex is the value of the (rows - 1)th row of the first occurence
				# firstoccurenceat is the index of the first occurence in the list PList
				# occurence is the total number of occurence of the pattern
				# gap is the gap between the consecutive occurences of the pattern

				if op(1..rows - 2, POS) = op(1..rows - 2, POSAdd) and 
						type(POS[rows], integer) and
						POSAdd[rows - 1] - POS[rows - 1] > lastgap then
					
					firstindex := POS[rows - 1]:
					firstoccurenceat := i;
					occurence := occurence + 1:
					gap := POSAdd[rows - 1] - POS[rows - 1]:
					resultlist := {(MIN_OCCURENCE - 1) * x1 + x2 = POSAdd[rows], (MIN_OCCURENCE - occurence) * x1 + x2 = POS[rows]};
				fi:
			else

				# find the next positions that might follow the pattern found
				if op(1..rows - 2, POS) = op(1..rows - 2, POSAdd) and 
						firstindex - POS[rows - 1] = gap then
					firstindex := POS[rows - 1]:
					occurence := occurence + 1:
					resultlist := resultlist union {(MIN_OCCURENCE - occurence) * x1 + x2 = POS[rows]};

#print(` FindPattern increasing occurence: `, ` POS = `, POS, ` firstindex = `, firstindex, ` firstoccurenceat = `, firstoccurenceat, ` occurence = `, occurence, ` gap = `, gap, resultlist);
					if occurence >= MIN_OCCURENCE then
						break;
					fi:
				fi:
			fi:
			
			# the pattern is invalid since we are unable to find the next desired position
			# loop back to the previous spot
#print(` FindPattern reset: `, ` POS = `, POS, ` firstindex = `, firstindex, ` firstoccurenceat = `, firstoccurenceat, ` occurence = `, occurence, ` gap = `, gap, resultlist);
			if i = 1 and gap <> 0 and firstoccurenceat <> 1 then
				gap := 0;
				occurence := 1;
				i := firstoccurenceat;
			fi:
		od:

		# remember the gap, so we know where to start next time around
		if gap <> 0 then 
			lastgap := gap;
		fi:
		
		# Got a possible pattern, verify it	
		if gap <> 0 and occurence >= MIN_OCCURENCE then
			sol := solve(resultlist, {x1, x2});
			
			if sol <> NULL and subs(sol, x1) >= 0 and subs(sol, x2) >= 0  then
#print(` FindPattern found: `, ` POS = `, POS, ` firstindex = `, firstindex, 
#	` firstoccurenceat = `, firstoccurenceat, ` occurence = `, occurence, 
#	` gap = `, gap, resultlist);
	
				patternlist := [op(1..rows - 2, POSAdd), 
					firstindex + gap * a[rows - 1], subs(sol, x1) * a[rows - 1] + subs(sol, x2)];

				# we found a solution, but it must not conflict with the existing nim values
				# i.e. it cannot coincide with "forbidden values" and 
				# it must agree with the existing patterns
				POS := patternlist;
					
#print(` find enough times: `, ` POS = `, POS, ` firstindex = `, firstindex, 
#	` firstoccurenceat = `, firstoccurenceat, ` occurence = `, occurence, ` gap = `, gap);

#if POSAdd = [1,3,44,50] and gap = 12then
#print(` Find Pattern in [1,3,44,50] PList = `, PList, ` NList = `, NList, gap, sol, POS, ElementsIntersect(POS, NList, a));
#fi:

				if not ElementsIntersect(POS, NList, a) 
						and not ElementsDisagree(POS, PList, a) then
PrintResult(cat(` Found SYMBOLIC position: `, POS, ` generated by `,
	seq( subs({a[rows- 1] = i, a = i}, POS), i = 0..MIN_OCCURENCE - 1)));

if [POS[1], POS[2]] = [4, 0] and POS[3] = 32*a[3]+32 then
print(`Why is this? `, POS, PList, NList);
fi:
					redudantlist := FindRedundantPositions(POS, PList, a);
					
					RETURN( POS, sort( convert(convert(PList, set) minus redudantlist[1]
						 union { redudantlist[2] }, list), SortChomp) );
				fi:
#print(` ElementsIntersect = `, ElementsIntersect(POS, NList, a), ` ElementDisagree = `, ElementsDisagree(POS, PList, a));
			else
#print(` FindPattern reset: `, ` POS = `, POS, ` firstindex = `, firstindex, ` firstoccurenceat = `, firstoccurenceat, ` occurence = `, occurence, ` gap = `, gap, resultlist);
				occurence := 1;
				firstoccurenceat := firstoccurenceat - 1;
			fi:
		fi:
		
#print(` FindPattern nothing found: `, ` POSAdd = `, POSAdd, ` POS = `, POS, ` firstindex = `, firstindex, ` firstoccurenceat = `, firstoccurenceat, ` occurence = `, occurence, ` gap = `, gap, ` lastgap = `, lastgap, resultlist);

		# we exaulsted all possibilities, so move on
		if lastgap >= len / (MIN_OCCURENCE - 1) or gap = 0 or len < MIN_OCCURENCE then
#print(` FindPattern nothing found: `, ` POS = `, POS, ` firstindex = `, firstindex, ` firstoccurenceat = `, firstoccurenceat, ` occurence = `, occurence, ` gap = `, gap, resultlist);
			break;
		fi:
	od:
#print(` FindPattern found patterns: `, patternlist);
				
	if bAddPOSIfNoPatternFound then
PrintResult(cat(` Found SINGLE position: `, POSAdd)):
		POSAdd, sort( [op(PList), POSAdd], SortChomp);
	else
		POSAdd, PList;
	fi:
end:

##############################################################################
# Find all the positions in PList that POS as a pattern already represents
# PList is assumed ordered
#
# return two elements
# the first is a list of redundant positions, solid or symbolic
# the sencond is the new formula found
#
# POS is a pattern
# PList is an ordered list of nim values
# a is the symbolic variable
##############################################################################

FindRedundantPositions := proc(POS, PList, a)
	local i, j, len, POS1, POS2, rows, bSolid, sol, returnset, newpattern;
	
	rows := nops(POS) - 1;
	len := nops(PList);
	bSolid := true;
	returnset := {};
	newpattern := POS;
	
	for i from len by -1 to 1 do
		# Search the list backward because we alread assumed PList is ordered
		
		POS1 := PList[i];
		
		if op(1..rows - NON_FIXED_ROWS, newpattern) = op(1..rows - NON_FIXED_ROWS, POS1) then 
			if degree(POS1[rows - NON_FIXED_ROWS + 1]) > 0 then
				# if POS1 is symbolic, then the two formulae shall match
				
				POS2 := subs(a = aa, newpattern);
				sol := solve({POS1[rows - NON_FIXED_ROWS + 1] - POS2[rows - NON_FIXED_ROWS + 1],
					POS1[rows + 1] - POS2[rows + 1]}, 
					{a[rows - NON_FIXED_ROWS + 1], aa[rows - NON_FIXED_ROWS + 1]});

				if sol = NULL then
					next;	
				elif (subs(sol, a[rows - NON_FIXED_ROWS + 1]) = a[rows - NON_FIXED_ROWS + 1] and 
						degree(subs(sol, aa[rows - NON_FIXED_ROWS + 1])) > 0 and 
						type(coeff(subs(sol, aa[rows - NON_FIXED_ROWS + 1]), a[rows - NON_FIXED_ROWS + 1], 0), integer)) or
					(subs(sol, aa[rows - NON_FIXED_ROWS + 1]) = aa[rows - NON_FIXED_ROWS + 1] and 
						degree(subs(sol, a[rows - NON_FIXED_ROWS + 1])) > 0 and
						type(coeff(subs(sol, a[rows - NON_FIXED_ROWS + 1]), aa[rows - NON_FIXED_ROWS + 1], 0), integer)) then
					
					if coeff(POS1[rows - NON_FIXED_ROWS + 1], a[rows - NON_FIXED_ROWS + 1], 0)
							= coeff(POS[rows - NON_FIXED_ROWS + 1], a[rows - NON_FIXED_ROWS + 1], 0) then
						if coeff(POS1[rows - NON_FIXED_ROWS + 1], a[rows - NON_FIXED_ROWS + 1], 1)
								> coeff(newpattern[rows - NON_FIXED_ROWS + 1], a[rows - NON_FIXED_ROWS + 1], 1) then
							POS2 := POS1;
						else
							POS2 := newpattern;
							newpattern := POS1;
						fi:	
					elif coeff(POS1[rows - NON_FIXED_ROWS + 1], a[rows - NON_FIXED_ROWS + 1], 0)
							> coeff(newpattern[rows - NON_FIXED_ROWS + 1], a[rows - NON_FIXED_ROWS + 1], 0) then
						POS2 := POS1;
					else
						POS2 := newpattern;
						newpattern := POS1;
					fi:	

					returnset := returnset union {POS2};
				fi:			
			elif bSolid then
				# if POS1 is solid, then
				POS2 := (POS1[rows - NON_FIXED_ROWS + 1] - coeff(newpattern[rows - NON_FIXED_ROWS + 1], a[rows - NON_FIXED_ROWS + 1], 0) ) 
						/ coeff(newpattern[rows - NON_FIXED_ROWS + 1], a[rows - NON_FIXED_ROWS + 1], 1);
						
				if type( POS2, integer ) then
					if POS2 >= 0 then
						if subs(a[rows - NON_FIXED_ROWS + 1] = POS2, POS1[rows + 1] - newpattern[rows + 1]) = 0 then
							# we found a larger position, erase it
							returnset := returnset union {POS1};
						fi:
					else
						if subs(a[rows - NON_FIXED_ROWS + 1] = POS2, POS1[rows + 1] - newpattern[rows + 1]) = 0 then
							# a smaller position, erase it
							returnset := returnset union {POS1, newpattern};
							newpattern := subs(a[rows - NON_FIXED_ROWS + 1] = a[rows - NON_FIXED_ROWS + 1] + POS2, newpattern);
						else
							# the position fits the formula, but the nim value does not
							# stop here, no more solid positions to be deleted
							bSolid := false;
						fi:
					fi:
				fi:
			fi:
		fi:
	od:

#if resultset <> {} then
#PrintResult(cat(` Redudant result deleted: `, returnset)); #, POS, PList)):
#fi:
	returnset, newpattern;
end:

##############################################################################
# Check that a nim value agrees with the values in LIST, useful exspecially
# if POS is symbolic
# 
# returns true if POS does not agree with some value if LIST
#
# POS is a Chomp position with its nim value
# LIST is the list of nim values to be checked against
# a is the symbolic variable
##############################################################################

ElementsDisagree := proc(POS, LIST, a)
	local i, j, len, POS1, POS2, aa, rows, eq, res, 
		var, bNotAgree, coef, coefa, coefaa, const;
	
	len := nops(LIST);
	rows := nops(POS) - 1:
	POS1 := subs(a = aa, POS);
		
#print(`ElementsDisagree: POS1 = `, POS1, `LIST = `, LIST):
	bNotAgree := false;
	
	for i from 1 to len do
		POS2 := LIST[i];
		
		bNotAgree := true:
#print(`ElementsDisagree: POS1 = `, POS1, `POS2 = `, POS2, `res = `, res):
		
		# The top rows must be the same
		for j from 1 to rows - NON_FIXED_ROWS do
			var := POS1[j] - POS2[j];
			if type(var, numeric) and var <> 0 then
				bNotAgree := false:
				break;
			fi:
		od:
		
		if not bNotAgree then
			next;
		fi:
		
		#!!!! only works if NON_FIXED_ROWS = 1 !!!!
		for j from rows - NON_FIXED_ROWS + 1 to rows do
			if EquationsDisagree(POS1[j] - POS2[j], POS1[rows+1] - POS2[rows+1], [a[j], aa[j]])
					or EquationsDisagree(POS1[j] - POS2[j], POS1[rows+1] - POS2[rows+1], [aa[j], a[j]])
					or EquationsDisagree(POS1[rows+1] - POS2[rows+1], POS1[j] - POS2[j], [a[j], aa[j]])
					or EquationsDisagree(POS1[rows+1] - POS2[rows+1], POS1[j] - POS2[j], [aa[j], a[j]]) then
				RETURN(true);
			fi:
		od:
	od:
	
	false;
end:

##############################################################################
# Find out if
##############################################################################

EquationsDisagree := proc(eq1, eq2, var)
	local i, j, a, aa, sol, temp, temp1;
	
	a := var[1];
	aa := var[2];
	
	if type(eq1, numeric) then
		# if eq1 is numeric but not 0, the two formulae never disagree with each other
		# but if eq1 is 0, then eq2 must be 0
		if eq1 <> 0 then
			RETURN(false):
		elif type(eq2, numeric) then
			if eq2 <> 0 then
				RETURN(true);
			else
				RETURN(false);
			fi:
		else
			RETURN(true);
		fi:
	else
		# if eq1 is not numeric, then if we can solve eq2 by the variables,
		# then there are other values for the variables to make the formulae disagree
		# so eq2 MUST be 0 in this case
		sol := solve({eq1}, {a});
		if sol = NULL then
			RETURN(false);
		else
			temp := subs(sol, a);
			
			if type(temp, numeric) then
				if not type(temp, integer) or temp < 0 then
					RETURN(false);	
				elif subs(sol, eq2) <> 0 then
					RETURN(true);
				else
					RETURN(false);
				fi:					
			else
				if not type(coeff(temp, aa, 0) * denom(coeff(temp, aa, 1)), integer) then
					RETURN(false);
				else
					sol := subs(sol, eq2);
					
					if type(sol, numeric) then
						if sol <> 0 then
							RETURN(true):
						else
							RETURN(false);
						fi:
					else
						RETURN(true);
					fi:
				fi:
			fi:
		fi:
	fi:
end:

##############################################################################
# Check if a nim value is part of the list LIST
# both POS and LIST might be symbolic, in that case
# we check if POS and LIST have an intersection point
#
# return true if POS agrees with some LIST at some point
#
# POS is a position with its nim value
# LIST is a list of "forbidden nim values" to be checked against
# a is a list of symbolic variables
##############################################################################

ElementsIntersect := proc(POS, LIST, a)
	local i, j, len, POS1, POS2, aa, rows, eq, res, 
		var, bIntersect, coef, coefa, coefaa, const;
	
	len := nops(LIST);
	rows := nops(POS) - 1:
	POS1 := subs(a = aa, POS):
	
#print(`ElementsIntersect: POS1 = `, POS1, `LIST = `, LIST):
	bIntersect := false;
	
	for i from 1 to len do
		POS2 := LIST[i];
		
		bIntersect := true:
#print(`ElementsIntersect: POS1 = `, POS1, `POS2 = `, POS2, `res = `, res):
			
		# The top rows must be the same
		for j from 1 to rows - NON_FIXED_ROWS do
			var := POS1[j] - POS2[j];
			if type(var, numeric) and var <> 0 then
				bIntersect := false:
				break;
			fi:
		od:
		
		if not bIntersect then
			next;
		fi:
		
		#!!!! only works if NON_FIXED_ROWS = 1 !!!!
		for j from rows - NON_FIXED_ROWS + 1 to rows do
			var := POS1[j] - POS2[j];
#print(`ElementsIntersect: var = `, var, `j = `, j, POS1, POS2);

			if type(var, numeric) then
				if var <> 0 then
					bIntersect := false:
					break;
				fi:
			else
				# if the last NON_FIXED_ROWS rows and the nim values agree, error
				var := solve({var, POS1[rows+1] - POS2[rows+1]}, {aa[j], a[j]});
#print(`ElementsIntersect: var = `, var, `j = `, j);
				if var = NULL then 
					bIntersect := false:
					break;
				fi:
				
				coefaa := subs(var, aa[j]);
				coefa  := subs(var, a[j]);
				if coeff(POS1[j], b, 1) <> 0 or coeff(POS2[j], b, 1) <> 0 then
					var := solve({coeffaa, coefa}, {b});
					var := subs(var, b);
					if type(var, integer) and var >= 0 then
						bIntersect := true;
					else
						bIntersect := false;
						break;
					fi:
				else
					# if there is a positive solution, then the two positions intersect
					if ((type(coefaa, integer) and coefaa >= 0) and (type(coefa, integer) and coefa >= 0)) or
							(coefaa = aa[j] and (not type(coefa, integer) and type(coeff(coefa, aa[j], 0), integer) and coeff(coefa, aa[j], 0) + coeff(coefa, aa[j], 1) > 0)) or
							(coefa = a[j] and (not type(coefaa, integer) and type(coeff(coefaa, a[j], 0), integer) and coeff(coefaa, a[j], 0) + coeff(coefaa, a[j], 1) > 0)) then
						bIntersect := true;
					else
						bIntersect := false;
						break;
					fi:
				fi;
#print(`ElementsIntersect: POS1 = `, POS1, `POS2 = `, POS2, `coefaa = `, coefaa, `coefa = `, coefa, `j = `, j);
			fi:
		od:
#print(`ElementsIntersect: POS1 = `, POS1, `POS2 = `, POS2, `res = `, res, `bIntersect = `, bIntersect):
		
		if bIntersect then
#if POS = [1,3,20+12*a[3],26+12*a[3]] then
#print(`ElementsIntersect: POS1 = `, POS1, `POS2 = `, POS2, `bIntersect = `, bIntersect):
#fi:
			break;
		fi:
	od:
	
	bIntersect:
end:

##############################################################################
# Build a "generating function" for calculation based on all the nim values
# and "forbidden nim values" derived from T for all positions 
# with fixed top rows FixedTop
#
# returns two elements
# the first one is the generating funciton
# the second one is all the "forbidden nim values" with the FixedTop
#
# T is the list of all nim values
# FixedTop is the value of the top rows of the position that we are calculating
# max_cols is the maximum cols allowed for the top rows
# a is a list of symbolic variables
##############################################################################

BuildGeneratingFunctionForFixedTop := proc(T, FixedTop, max_cols, a)
	local i, j, k, l, ll, len, len1, POS, POS1, POS2, listT, LIST, lenLIST, lenlistT, 
		returnset, returnset2, set1, set2, rows, firstrow, f, f1, convf, NList;
	
	rows := nops(FixedTop) + NON_FIXED_ROWS;
	firstrow := FixedTop[1];
	NList := [];
	
#print(`BuildGeneratingFunctionForFixedTop: T = `, T, `  FixedTop = `, FixedTop, `  a = `, a);
	returnset := {};
	returnset2 := {};
	
	f := 0;
#print(`BuildGeneratingFunctionForFixedTop: f = `, f);

	if nops(T) >= rows then
		len := nops(T[rows]):
		if len >= firstrow then
			len := firstrow - 1:
		fi:
#print(`BuildGeneratingFunctionForFixedTop: len = `, len, `firstrow = `, firstrow):
		
		# Get "forbidden nim values" from all known nim values with the same number of rows
		
		# check "forbidden nim values" with the first row less than this one
		for i from 1 to len do:
			LIST := T[rows][i];
			lenLIST := nops(LIST);
			for j from 1 to lenLIST do
				POS := LIST[j]:
				POS1 := [op(1..nops(POS) - 1, POS)];
				set1 := GrowPositionFromTopByMaxValue(POS1, firstrow - i, max_cols, a):

				set2 := set1[2];
				set1 := set1[1];
				
				# calculate all the positions grown from the top, but with limited expansion
				# i.e. only the top row has been changed
				len1 := nops(set1);
				for k from 1 to len1 do
					if [op(1..rows - NON_FIXED_ROWS, set1[k])] = FixedTop then
						convf := ConvertPositionToGeneratingFunctionWithFixedTop(set1[k], POS, a):
						NList := [op(NList), op(convf[2])];
						f := f + convf[1];
					fi:
				od:
				
				# calculate all the positions grown from the top, but with unlimited expansion
				# i.e. the top row used to be the same as the bottom row
				# so, no furthur positions can have the nim value any more 
				len1 := nops(set2);
				if len1 <> 0 then
					for k from 1 to len1 do
						if [op(1..rows - NON_FIXED_ROWS, set2[k])] = FixedTop then
							POS2 := POS[nops(POS)];
							
							for l from 1 to NON_FIXED_ROWS do
								POS2 := subs({a[nops(FixedTop) + l]=0}, POS2):
							od:
							
							POS2 := [op(1..nops(FixedTop), POS), seq(0, l=1..NON_FIXED_ROWS), POS2];

							convf := ConvertPositionToGeneratingFunctionWithFixedTop
								(set2[k], POS2, a):
							NList := [op(NList), op(convf[2])];
							f := f + convf[1];
						fi:
					od:
				fi:
			od:
		od:
		
		# check "forbidden nim values" with the first row the same as this one
		if nops(T[rows]) >= firstrow then
			LIST := T[rows][firstrow];
			lenLIST := nops(LIST);
#print(`BuildGeneratingFunctionForFixedTop: listT = `, LIST):
			
			for j from 1 to lenLIST do
				POS := LIST[j];
				POS1 := [op(1..nops(POS) - 1, POS)];
				
				if [op(1..rows - NON_FIXED_ROWS, POS1)] = FixedTop then
					# if the top is the same as the FixedTop, do not expand it
					# it is because the positions with the values of their last row as 0
					convf := ConvertPositionToGeneratingFunctionWithFixedTop(POS1, POS, a):
					NList := [op(NList), op(convf[2])];
					f := f + convf[1];
				else
					set1 := GrowPositionFromBottomByMaxValue(POS1, max_cols, a):
					set2 := set1[2];
					set1 := set1[1];
#print(`BuildGeneratingFunctionForFixedTop: listT = `, listT, POS, POS1, set1, set2):
					
					# calculate all the positions grown from the bottom, if the last row is not 0
					len1 := nops(set1);
					for k from 1 to len1 do
						if [op(1..rows - NON_FIXED_ROWS, set1[k])] = FixedTop then
							convf := ConvertPositionToGeneratingFunctionWithFixedTop(set1[k], POS, a):
							NList := [op(NList), op(convf[2])];
							f := f + convf[1];
						fi:
					od:

					# calculate all the positions grown from the bottom, if the last row is 0
					len1 := nops(set2);
					if len1 <> 0 then
						for k from 1 to len1 do
							if [op(1..rows - NON_FIXED_ROWS, set2[k])] = FixedTop then
								POS2 := POS[nops(POS)];
								
								for l from 1 to NON_FIXED_ROWS do
									POS2 := subs({a[nops(FixedTop) + l]=0}, POS2):
								od:
								
								POS2 := [op(1..nops(FixedTop), POS), seq(0, l=1..NON_FIXED_ROWS), POS2];

								convf := ConvertPositionToGeneratingFunctionWithFixedTop
									(set2[k], POS2, a):
								NList := [op(NList), op(convf[2])];
								f := f + convf[1];
							fi:
						od:
					fi:
				fi:
			od:
		fi:
	fi:
#print(`BuildGeneratingFunctionForFixedTop: firstrow = `, firstrow, ` returnset = `, returnset);

	# Get "forbidden nim values" from all positions with less rows

	for i from 1 to rows - 1 do
		listT := T[i]:
		len := nops(listT):
#print(`BuildGeneratingFunctionForFixedTop: listT = `, listT):

		for j from 1 to len do
			len1 := nops(listT[j]);

			for k from 1 to len1 do
				POS := listT[j][k];
				POS1 := [op(1..nops(POS) - 1, POS)];
#print(`BuildGeneratingFunctionForFixedTop: listT = `, listT, POS, POS1):

				set1 := ExpandPositionWithNewRow(POS1, rows, firstrow, max_cols, a):

				# calculate all the positions grown from positions with less rows
				for l from 1 to nops(set1) do
					if	[op(1..rows - NON_FIXED_ROWS, set1[l])] = FixedTop then
						convf := ConvertPositionToGeneratingFunctionWithFixedTop(set1[l], POS, a):
						NList := [op(NList), op(convf[2])];
						f := f + convf[1];
					fi:
				od:
			od:
#print(`BuildGeneratingFunctionForFixedTop: listT = `, listT, `j = `, j, `LIST = `, LIST):

		od:	
	od:
#print(`BuildGeneratingFunctionForFixedTop: firstrow = `, firstrow, ` returnset = `, returnset);

	f1 := 1;
	for i from 1 to NON_FIXED_ROWS do
		f1 := f1 / (1 - Symb_Var[nops(FixedTop) + i]);
	od:
	f1 := f1 / (1 - Symb_Nim);

	f1 - f, NList;	
end:

#####################################################################################
# This function is to return the set of all the positions grown from POS 
# by adding new rows. 
# So any positions grown from POS cannot have the same nim value as POS does
#
# return a set of positions expanded from POS
#
# POS is the original position
# rows is the number of rows the new positions are to have
# firstrow is the value of the first row of the new position
# a is a list of symbolic variables
#####################################################################################

ExpandPositionWithNewRow := proc(POS, rows, firstrow, max_cols, a)
	local i, len, NEWPOS;
	
	len := nops(POS):

	# reorganize a, i.e. a[i] --> a[i + rows - len] in POS
	NEWPOS := [seq(0, i=1..rows-len), 
		seq(coeff(POS[i], a[i], 0) + coeff(POS[i], a[i], 1) * a[rows - len + i], i = 1..len)];

	GrowPositionFromTopByMaxValue(NEWPOS, firstrow, max_cols, a)[1]:		
end:


##############################################################################
# Convert a position into generating function, 
# i.e. the monomial representation of the Chomp position
#
# returns two elements
# the first one is the generating function
# the second one is the "forbidden nim values" generated from POS
#
# POS is a Chomp position
# ParentPOS is the position that POS is generated from
# a is the symbolic variable
##############################################################################

ConvertPositionToGeneratingFunctionWithFixedTop := proc(POS, ParentPOS, a)
	local i, len, lenparent, nimvalue, f, f1, denomf, diffa, NList, POS1;

#print(`ConvertPositionToGeneratingFunctionWithFixedTop: POS = `, POS, `ParentPOS = `, ParentPOS, `f = `, f):

	f := 1:
	denomf := 1:
	len := nops(POS);
	lenparent := nops(ParentPOS) - 1;
	nimvalue := ParentPOS[lenparent + 1];
	
	POS1 := [seq(0, i=1..len - lenparent)];
	for i from 1 to lenparent do
		POS1 := [op(POS1), subs(a[i] = a[i + len - lenparent], ParentPOS[i])];
		nimvalue := subs(a[i] = a[i + len - lenparent], nimvalue);
	od:
	
	# we need to know how much have been drawn from the symbolic variable
	# to make the expansion
	diffa[len - NON_FIXED_ROWS] := 0;
	for i from 1 to len - NON_FIXED_ROWS do
		diffa[len - NON_FIXED_ROWS] := diffa[len - NON_FIXED_ROWS] + POS[i] - POS1[i];
	od;
	for i from len - NON_FIXED_ROWS + 1 to len do
		diffa[i] := diffa[i-1] + POS[i] - POS1[i];
	od;
	
	# adjust the values by the coefficients of the symbolic variables
	# because for every number added to the POS ????
	for i from len - NON_FIXED_ROWS + 1 to len do
		if coeff(POS1[i], a[i], 1) = 0 then
			if diffa[i] <> 0 and POS1[i] <> 0 then
				ERROR(`The result does not look right for conversion: from `, POS1, ` to `, POS);
			fi:
		else
			diffa[i] := diffa[i] / coeff(POS1[i], a[i], 1);
			if not type(diffa[i], integer) then
				ERROR(`The result does not look right for conversion: from `, POS1, ` to `, POS);
			else
				diffa[i] := diffa[i] * coeff(nimvalue, a[i], 1);
			fi;
		fi:
	od;	
	
	if [op(1..len - NON_FIXED_ROWS, POS)] <> [op(1..lenparent - NON_FIXED_ROWS, ParentPOS)] then
		# if the positions is generated from positions with lesser rows,
		# or if they are of the same number of rows, but not with the same top rows
		# we only have to consider the position itself
		
		# the nim value cannot be any of the ones from "smaller" positions
		# and cannot be any of the ones with lesser last row

		for i from len - NON_FIXED_ROWS + 1 to len do
			f := f * Symb_Var[i] ^ coeff(POS[i], a[i], 0):
			denomf := denomf * Symb_Var[i] ^ coeff(POS[i], a[i], 1):
		od:
		
		if not type(ParentPOS[lenparent+1], integer) then
			f := f * Symb_Nim ^ (coeff(nimvalue, a[len], 0) + diffa[len]):
		else
			f := f * Symb_Nim ^ (coeff(nimvalue, a[len], 0)):
		fi:
		denomf := denomf * Symb_Nim ^ coeff(nimvalue, a[len], 1):
		
		if not type(denomf, integer) then
			f := f / (1 - denomf);
		fi:
		
		if len = lenparent and ParentPOS[lenparent] = 0 then
			NList := [[op(1..len - NON_FIXED_ROWS, POS), b, nimvalue]];
		else
			if type(ParentPOS[lenparent], integer) and diffa[len] <> 0 then
				ERROR(`what happened in ConvertPositionToGeneratingFunctionWithFixedTop????`)
			fi:
			NList := [[op(1..len - NON_FIXED_ROWS, POS), POS[len], nimvalue + diffa[len]]];
		fi;	
	else
		# the nim value cannot be the same as any of the ones from "smaller" positions
		# once a position has a nim value, all the bigger number shall be elminated too
		
		f1 := 0;
		for i from len - NON_FIXED_ROWS + 1 to len do
			f := f * Symb_Var[i] ^ coeff(POS[i], a[i], 0):
			denomf := denomf * Symb_Var[i] ^ coeff(POS[i], a[i], 1):
			f1 := f1 + 1 / (1 - Symb_Var[i]);
		od:
		f := f * Symb_Nim ^ coeff(nimvalue, a[len], 0):
		denomf := denomf * Symb_Nim ^ coeff(nimvalue, a[len], 1):
		f1 := f1 + 1 / (1- Symb_Nim);
		
		f := f * f1;
		if not type(denomf, integer) then
			f := f / (1 - denomf);
		fi:
		NList := [[op(1..len - NON_FIXED_ROWS, ParentPOS), ParentPOS[len] + b, nimvalue]];
	fi:

#if [POS[1], POS[2]] = [1, 1] then
#print(`ConvertPositionToGeneratingFunctionWithFixedTop: POS = `, POS, `ParentPOS = `, ParentPOS, `f = `, f, ` NList = `, NList):
#fi:

	f, NList:
end:

#####################################################################################
# POS is a Chomp position, grow POS to POS1 so that POS1[1] = POS[1] + c
# and POS can be deduced from POS1 by one move
# which means if POS is nim value, then the result is the set of 
# "forbidden nim values" grown by POS with extra pieces at the top
# NOTE: POS itself might be in the result
#
# return two elements
# the first is the expansion from POS, if the second row is not 0
# the second is the expansion from POS, if only the first row is not 0
# they are treated differently because the generating function 
# are not the same in the two occasions
#
# POS is a Chomp position
# c is a positive integer
# max_cols is the maximum cols allowed for the top rows
# a is a list of symbolic variables
#####################################################################################

GrowPositionFromTopByMaxValue := proc(POS, c, max_cols, a) 
    local i, j, ret:
#print(`GrowPositionFromTopByMaxValue: POS = `, POS, `c = `, c):

	if nops(POS) = 1 then
		RETURN({[POS[1] + c]}, {});
    elif type(POS[2], integer) then
        if c <= POS[2] then			# only change/grow the first row
            RETURN({[POS[1] + c, POS[2] - c, seq(POS[i], i = 3..nops(POS))]}, {}):

        elif POS[2] = 0 then		# we have to change a block of 0s
            for i from 2 to nops(POS) do	# find the size of the 0 block
                if not(type(POS[i], integer) and POS[i] = 0) then
					break;
                fi:
            od:

			# now j is the last row of the 0 block
			j := i - 1;
            if j < nops(POS) then
				if coeff(POS[j+1], a[j+1], 0) < c then
				    if coeff(POS[j+1], a[j+1], 1) > 0 then	
				    
						# if the row has a symbol, extend beyond it
						# be careful the coeff of POS[j + 1] can be > 1
				        ret := ExpandBlockByMaxValue([POS[1] + c], max_cols, 
							[coeff(POS[j+1], a[j+1], 1) * a[j+1] 
							+ modp(coeff(POS[j+1], a[j+1], 0) - c, coeff(POS[j+1], a[j+1], 1)), 
							seq(POS[i], i = j+2..nops(POS))], j - 1, a):
						
						# !!!! works only when NON_FIXED_ROWS = 1 !!!!
						if coeff(POS[j+1], a[j+1], 0) = 0 then
							# even if the row can be symbolic, it still can be treated as 0
							RETURN( ret, GrowPositionFromTopByMaxValue([op(1..j, POS), 0], c, max_cols, a)[2]);
						else
							RETURN( ret, {});
						fi:
				    else									# cannot grow the position
				        RETURN( {}, {} ):
				    fi:
				fi:

				# expand the block
				RETURN(ExpandBlockByMaxValue([POS[1] + c], max_cols, 
					[POS[j+1] - c, seq(POS[i], i = j+2..nops(POS))], j - 1, a), {}):
			else
				# all the rows are 0 except the first one, 
				# so the result is the second item
				RETURN( {}, ExpandBlockByMaxValue([POS[1] + c], max_cols, 
					[a[j]], j - 2, a) ):
			fi:
        else	# c > POS[2] and POS[2] <> 0, cannot grow the position
            RETURN( {}, {} ):
        fi:

    fi:

	# now POS[2] is symbolic, then POS has two rows
	# !!!! works only when NON_FIXED_ROWS = 1 !!!!
	if nops(POS) > 2 then
		ERROR(`Unable to handle input in GrowPositionFromTopByMaxValue: `, POS);
	else								# two rows
		if coeff(POS[2], a[2], 0) >= c then		# same as the instance above
			RETURN({[POS[1] + c, POS[2] - c]}, {}):
		elif coeff(POS[2], a[2], 0) < c and coeff(POS[2], a[2], 1) > 0 then	
			# expand beyond the symbol
		    ret := {[POS[1] + c, coeff(POS[2], a[2], 1) * a[2] 
				+ modp(coeff(POS[2], a[2], 0) - c, coeff(POS[2], a[2], 1))]}:
			
			if coeff(POS[2], a[2], 0) = 0 then
				# if only the first row is not 0
				RETURN(ret, {[POS[1] + c, a[2]]});
			else
				RETURN(ret, {});
			fi:
		else 
			ERROR(`UNEXPECTED RESULT: `, POS);
		fi:
	fi:
end:

#####################################################################################
# Grow a position POS to POS1 so that POS[1] = POS1[1] 
# and POS can be deduced from POS1 by one cut
# which means if POS is nim value, then the result is the set of 
# "forbidden nim values" grown by POS with the same number of rows 
# and the same first row
# NOTE: POS itself might be in the result
#
# return two elements
# the first is the expansion from POS, if the last rows are not 0
# the second is the expansion from POS, if the last rows are 0
# they are treated differently because the generating function 
# are not the same in the two occasions
#
# POS is the original position
# max_cols is the maximum cols allowed for the top rows
# a is a list of symbolic variables
#####################################################################################

GrowPositionFromBottomByMaxValue := proc(POS, max_cols, a)
	local i, j, k, jj, value, para, len, list3, list4, list1, list2:
    
    list3 := {}: # the regular return value
    list4 := {}; # the return value for last row = 0
    len := nops(POS):
#print(`GrowPositionFromBottomByMaxValue: POS = `, POS):
    
    for i from 2 to len - 1 do	#expand at each rows separately
        value := coeff(POS[i + 1], a[i + 1], 0):
        para := coeff(POS[i + 1], a[i + 1], 1):
        
        if value = 0 and para = 0 then
			# If the next row is 0, be careful. Find out how big the block is
			for j from i + 1 to nops(POS) do	# find out how many rows are 0
				if not (type(POS[j], integer) and POS[j] = 0) then
					break;
				fi:
			od:
			k := j - 1;
			
			list2 := [seq(POS[j], j = 1..i-1)]:
			
			# Expand the rows, including the current row treated as 0
			list1 := convert(ExpandBlockByMaxValue(list2, max_cols, 
				[seq(POS[j], j = k + 1..len)], k - i + 1, a), list):
			
			# Put the value of the current row back into the list
			list1 := {seq( [ op(list2), 
				POS[i] + op(i, list1[jj]), 
				seq(list1[jj][j], j=i+1..len)], jj = 1..nops(list1)) }:
#print(`GrowPositionFromBottomByMaxValue: list1 = `, list1, ` POS = `, POS);

			if k = nops(POS) then
				# the last row is 0
				list4 := list1;
			else
				list3 := list3 union list1:
			fi:
						
			i := k: # now we can skip the rows we have calculated
        else
			# otherwise, just expand the current row corresponding to the value of the next row
			list3 := list3 union { seq( [seq(POS[j], j = 1..i-1), POS[i] + k, POS[i + 1] - k,  
										seq(POS[j], j = i + 2..len)], k = 1..value ) }:
#print(`GrowPositionFromBottomByMaxValue: POS = `, POS):

			if para > 0 then # in this case, i <> len - 1 because POS shall be a nim value
				# extend the row to max_cols so that the expansion can cover all possible occasions
			    list3 := list3 union { seq( [op(1..i-1, POS), 
					POS[i] + value + j, para * a[i+1] + modp(para - j, para),  
					op(i + 2..len, POS) ], j = 1..max_cols - value ) }:
			fi:
		fi:
    od:

	if coeff(POS[len], a[len], 1) <> 0 and coeff(POS[len], a[len], 0) = 0  then
		# if the last row is symbolic, then we treat it as 0 too
		list4 := GrowPositionFromBottomByMaxValue([op(1..len-1, POS), 0], max_cols, a);
		list4 := list4[2];
	fi:
	
	list3, list4;
end:


#####################################################################################
# Expand a group of rows with the same size 
#
# If the original position looks like [prefix[1],..., prefix[a], 0,..., 0, suffix[1],..., suffix[b]]
# the result shall look like [prefix[1],..., prefix[a], alpha,..., beta, 
#                                  suffix[1] - alpha - ... - beta, suffix[2],..., suffix[b]]
# so that the result can be deduced to the original position by one cut
# the value of alpha,...,beta are determined by the value of suffix[1]
# NOTE: it might return the original position too, but we don't care
#
# prefix is a list that is preceding the result
# max_value is the maximum value that the rows can grow to
# suffix is a list that follows the result
# rows is the size of the block
# a is a list of symbolic variables
#####################################################################################

ExpandBlockByMaxValue := proc(prefix, max_value, suffix, rows, a)
    local i, CoeSymb, value, sum_to_value, list1, list2, POS:

#print(`ExpandBlockByMaxValue prefix = `, prefix, `max_value = `, max_value, `suffix = `, suffix, `rows = `, rows); 

    value := max_value - sum(prefix[i], i = 1..nops(prefix));
    
    if value < 0 then
		RETURN( {} );
    elif nops(suffix) = 0 then # if this is the last block, no restiction whatsoever

		# just expand to the second to the last row of all possible values, 
		# but we must maintain 0 as the constant part of the value in the last row, 
		# because of ConvertPositionToGeneratingFunctionWithFixedTop
        CoeSymb := 1:
        value := value;
        sum_to_value := value;
		
        list1 := convert(RecursiveExpandBlockByMaxValue(prefix, value, sum_to_value,
				[value], rows-1, CoeSymb, a), list);

		value := nops(list1);
		sum_to_value := nops(prefix) + rows;
		
		list2 := [];
		for i from 1 to value do
			POS := list1[i];
			POS := [op(1..sum_to_value - 1, POS), a[sum_to_value] * coeff(POS[sum_to_value], a[sum_to_value], 1)];
			
			list2 := [op(list2), POS];
		od:

		RETURN( list2 );
				
    elif coeff(suffix[1], a[nops(prefix) + rows + 1], 1) > 0 then 
		# we can symbolically expand the block
        CoeSymb := coeff(suffix[1], a[nops(prefix) + rows + 1], 1);
        sum_to_value := coeff(suffix[1], a[nops(prefix) + rows + 1], 0);
        value := value; 
    else
		# just expand to the value of the next row
        CoeSymb := 0:
        value := min(value, suffix[1]):
        sum_to_value := suffix[1];
    fi:
  
    # first row of suffix is replaced in RecursiveExpandBlockByMaxValue
    RecursiveExpandBlockByMaxValue(prefix, value, sum_to_value,
		[seq(suffix[i], i=2..nops(suffix))], rows, CoeSymb, a):
end:

#####################################################################################
# Recursive function used by ExpandBlockByMaxValue
#
# prefix is from ExpandBlockByMaxValue
# max_value is the maximum that the sum of all the new rows can reach
# sum_to_value is the minimum that the sum of all the new rows shall reach
# suffix is [suffix[2],...,suffix[b]] from ExpandBlockByMaxValue
# bExtendToInfinity is true if suffix[1] is symbolic
# a is a list of symbolic variables
#####################################################################################

RecursiveExpandBlockByMaxValue := proc(prefix, max_value, sum_to_value, suffix, rows, CoeSymb, a)
    local i:
#print(`RecursiveExpandBlockByMaxValue -- `, prefix, max_value, sum_to_value, suffix, rows, CoeSymb);  

	# we assume sum_to_value always >= max_value if old suffix[1] is not symbolic
	if sum_to_value < 0 and CoeSymb = 0 then
		# if the no symbolic value, sum_to_value cannot be negative
		RETURN( {} ):
    elif rows = 0 then
		RETURN({[op(prefix), a[nops(prefix)+1] ^ CoeSymb, op(suffix)]});
    elif rows = 1 then	# build the single row 
		# reconstruct the old suffix[1], but we only extend the last two rows to infinity
		
        if CoeSymb > 0 then
			# suffix has one row, so the last one expanded rows to infinity
			# remember: sum_to_value can be negative, which is to compensated by the symbolic value
			# this is represented in the second item of the union
				
			RETURN( {seq([op(prefix), i, CoeSymb * a[nops(prefix) + 2] + sum_to_value - i, 
					op(suffix)], i = 0..min(sum_to_value, max_value))}
				union {seq([op(prefix), i, CoeSymb * a[nops(prefix) + 2] 
					+ modp( - i + sum_to_value, CoeSymb), op(suffix)], 
					i = max(sum_to_value + 1, 0)..max_value)} ):
        else
			if nops(suffix) = 0 then
				# if this is the second to last row, we can add up to sum_value
				RETURN( {seq([op(prefix), i, sum_to_value - i], i = 0..sum_to_value)} ):
			else
				# if not, only up to max_value
				RETURN( {seq([op(prefix), i, sum_to_value - i, 
					op(suffix)], i = 0..min(sum_to_value, max_value))} ):
			fi:
        fi:
    else	
		# add the value of the first row in the block to the prefix, and run the function again
        RETURN( { seq(op(RecursiveExpandBlockByMaxValue([op(prefix), i], 
				max_value - i, sum_to_value - i,
				suffix, rows - 1, CoeSymb, a)), i = 0..max_value) } ):
    fi:
end:

##############################################################################
# Sort Chomp positions. This function is used by sort Maple function
# sort starts from the first to the last row
# return true if the value of a row of POS1 < that of POS2
# and all the previous rows are the same
#
# POS1 is a Chomp position
# POS2 is a Chomp position
##############################################################################

SortChomp := proc(POS1, POS2)
	local i, j, len1, len2, coef1, coef2;
	
	len1 := nops(POS1):
	len2 := nops(POS2):
	
	if len1 = len2 then
		for i from 1 to len1 - 1 do
			coef1 := coeff(POS1[i], a[i], 0):
			coef2 := coeff(POS2[i], a[i], 0):
			
			if coef1 < coef2 then
				RETURN( true ):
			elif coef1 > coef2 then
				RETURN(false):
			fi:
		od:
		
		RETURN(true);
	elif len1 < len2 then
		RETURN(true):
	else
		RETURN(false):
	fi:
end:

#####################################################################################
# Find the nim value of a Chomp position
#
# returns an integer as the nim value
#
# POS is the Chomp position whose nim value is to be found
# a is the symbolic variable used in InitNim
#####################################################################################

FindPositionNimValue := proc(POS, a)
	local nimlist, len, i, j, POS1, lenPOS, nimvalue, temp;
	
	lenPOS := nops(POS);
	
	if nops(InitNim) < lenPOS or nops(InitNim[lenPOS]) < POS[1] then
		ERROR(cat(`Input too large for the result set. Try any one with less rows or columns`, 
			nops(InitNim), lenPOS, nops(InitNim[lenPOS]), POS[1]));
	fi:
	
	nimlist := InitNim[lenPOS][POS[1]];
	len := nops(nimlist):
	
	for i from 1 to len do
		POS1 := nimlist[i];
		nimvalue := POS1[lenPOS + 1];
		
		if POS = [op(1..lenPOS, POS1)] then
			RETURN(nimvalue);
		elif degree(nimvalue) > 0 then
			# symblic solution
			
			if [op(1..lenPOS - NON_FIXED_ROWS, POS)] = [op(1..lenPOS - NON_FIXED_ROWS, POS1)] then
				for j from lenPOS - NON_FIXED_ROWS + 1 to lenPOS do
					if degree(POS1[j]) > 0 then
						temp := (POS[j] - coeff(POS1[j], a[j], 0)) / coeff(POS1[j], a[j], 1);
						
						if temp < 0 or not type(temp, integer) then
							break;
						fi:
						
						nimvalue := subs(a[j] = temp, nimvalue);
					elif POS1[j] <> POS[j] then
						break;
					fi:
				od:
				
				# the loop was finished, and the result is found
				if j = lenPOS + 1 then
					RETURN(nimvalue);
				fi:
			fi:
		fi:
	od:
	
	ERROR(`Input too large for the result set. Try any one with less rows or columns`);
end:

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

PrintResult := proc(msg)
	if PRINT_RESULT and printlevel > 1 then
		print(msg);
	fi;
end:

#####################################################################################
# Compare two results to find extra values in each one of them
#
# Nim1, Nim2 two values to be compared
#####################################################################################

CompareResult := proc(Nim1, Nim2)
	local i, j, len, len1, len2, top, top1, top2, diff1, diff2;
	
	len1 := nops(Nim1);
	len2 := nops(Nim2);
	len := min(len1, len2);
	
	for i from 1 to len do
		top1 := nops(Nim1[i]);
		top2 := nops(Nim2[i]);
		top := min(top1, top2);
		
		for j from 1 to top do
			diff1 := convert(Nim1[i][j], set) minus convert(Nim2[i][j], set);
			diff2 := convert(Nim2[i][j], set) minus convert(Nim1[i][j], set);
			
			if diff1 <> {} then
				print(cat(`Extra value in the FIRST: `, diff1));
			fi:
			if diff2 <> {} then
				print(cat(`Extra value in the SECOND: `, diff2));
			fi:
			if diff1 = {} and diff2 = {} then
				print(cat(`Results are the SAME for rows `, i, ` with top `, j));
			fi:	
		od:
		
		if top1 > top then
			print(cat(`Extra value in the FIRST: `, [op(top+1..top1, Nim1[i])]));
		fi:
		if top2 > top then
			print(cat(`Extra value in the SECOND: `, [op(top+1..top2, Nim2[i])]));
		fi:
	od:
	
	if len1 > len then
		print(cat(`Extra value in the FIRST: `, [op(len+1..len1, Nim1)]));
	fi:
	if len2 > len then
		print(cat(`Extra value in the SECOND: `, [op(len+1..len2, Nim2)]));
	fi:
end:

#####################################################################################
# Return a list of all monomials in a polynomial
#
# P is a polynomial
#####################################################################################

GetItemsFromPolynomial := proc(P)
	local i, list1;
	
	if type(P, `+`) then	# get all monomials as a list
		RETURN([op(P)]);
	else					# the polynomial is a monomial
		RETURN([P]):	
	fi:
end:

#####################################################################################
# Get all the positions with the desired nim value from InitNim
#
# return the list of the positions
#
# nim is the desired nim value
# a is a list of symbolic variables
#####################################################################################

GetPositionsWithNimValue := proc(nim, a)
	local i, j, k, list1, list2, listret, listret1, listret2, pos, rows, posnim, var;

	listret := [];
	
	for i from 1 to nops(InitNim) do
		list1 := InitNim[i];
		listret2 := [];
		
		for j from 1 to nops(list1) do
			list2 := list1[j];
			listret1 := [];
			
			for k from 1 to nops(list2) do
				pos := list2[k];
				
				rows := nops(pos) - 1;
				posnim := pos[rows + 1];
				
				if type(posnim, integer) then
					if posnim = nim then
						listret1 := [op(listret1), [op(1..rows, pos)]];
					fi:
				else
					var := solve(posnim = nim, a[rows]);
					if var <> NULL and var >=0 and type(var, integer) then
						listret1 := [op(listret1), subs(a[rows] = var, [op(1..rows, pos)])];
					fi:						
				fi:
			od:
			
			if listret1 <> [] then
				listret2 := [op(listret2), listret1];
			fi:
		od:
		
		if listret2 <> [] then
			listret := [op(listret), listret2];
		fi:
	od:

	listret;
end:

##############################################################################
# Get the next move so that the resulted position has the desired nim value 
#
# returns the next position if it can be found, NULL otherwise
#
# POS is a Chomp position
# nim is the desired nim value
# a is a list of symbolic variables
##############################################################################

FindNextMoveWithNim := proc(POS, nim, a)
	local i, j, k1, k2, k3, len1, len2, len3, rows, NList, NList1, FixedTop, 
		eq, res, var, POS2, max_cols, InitP, InitP1;

	InitP := GetPositionsWithNimValue(nim, a);
#print(InitP);
	
	rows := nops(POS):
	FixedTop := [op(1..rows - NON_FIXED_ROWS, POS)];

	max_cols := sum(POS[ii], ii = 1..rows - NON_FIXED_ROWS);
		
	for k1 from 1 to min(rows - 1, nops(InitP)) do
		len1 := nops(InitP[k1]);

		for k2 from 1 to len1 do
			len3 := nops(InitP[k1][k2]);
				
			for k3 from 1 to len3 do
				POS2 := InitP[k1][k2][k3];
				
				if POS = POS2 then 
					RETURN();
				fi:
				
				NList := ExpandPositionWithNewRow(POS2, rows, POS[1], max_cols, a);
#print(`FindWinningMove 1: POS2 = `, POS2, ` NList = `, NList);

				if ElementInList(POS, NList, a) then
#print(`FindWinningMove 1: POS2 = `, POS2, ` NList = `, NList);
					RETURN( POS2 );
				fi:
			od;
		od:
	od:
	
	if nops(InitP) >= rows then
		len1 := min(nops(InitP[rows]), POS[1] - 1);
		
		for k1 from 1 to len1 do
			len2 := nops(InitP[rows][k1]);

			for k2 from 1 to len2 do
				POS2 := InitP[rows][k1][k2];
				
				if POS = POS2 then 
					RETURN();
				fi:
				
				NList := GrowPositionFromTopByMaxValue(POS2, POS[1] - POS2[1], max_cols, a);
#print(`FindWinningMove 2: POS2 = `, POS2, ` NList = `, NList);
				NList := NList[1] union NList[2];
			
				if ElementInList(POS, NList, a) then
					RETURN( POS2 );
				fi:
			od:
		od:
		
		if nops(InitP[rows]) >= POS[1] then
			len1 := nops(InitP[rows][POS[1]]);
			
			for k1 from 1 to len1 do
				POS2 := InitP[rows][POS[1]][k1];
				
				if POS = POS2 then 
					RETURN();
				fi:
				
				NList := GrowPositionFromBottomByMaxValue(POS2, max_cols, a);
				NList := NList[1] union NList[2];
				# add extra element because it may not be in NList yet
				NList := NList union {[op(1..rows-1, POS2), POS2[rows] + a[rows]]};
#print(`FindWinningMove 3: POS2 = `, POS2, ` NList = `, NList);
			
				if ElementInList(POS, NList, a) then
					RETURN( POS2 );
				fi:
			od:
		else
PrintResult(cat(`FindWinningMove: Nothing found for `, POS, ` but the data available is incomplete`));
			RETURN();
		fi:
	fi:
	
PrintResult(cat(`FindWinningMove: Nothing found for `, POS));
	RETURN():
end:

##############################################################################
# Given a position, check if it is in a list
#
# return true if it is, false otherwise
#
# POS is a solid Chomp position
# NList is a list of solid and symbolic positions
# a is a list of symbolic variables
##############################################################################

ElementInList := proc(POS, NList, a)
	local rows, i, j, var, POS1;
#print(`ElementInList: POS = `, POS, ` NList = `, NList);
	
	rows := nops(POS);
	
	for i from 1 to nops(NList) do
		POS1 := NList[i];
		
		if [op(1..rows-1, POS1)] = [op(1..rows-1, POS)] then
			if type(POS1[rows], integer) then
#print(`ElementInList integer: POS = `, POS, ` NList = `, NList);
				if POS1[rows] = POS[rows] then
					RETURN(true);
				else
					RETURN(false);
				fi:				
			else
#print(`ElementInList not integer: POS = `, POS, ` NList = `, NList);
				var := solve(POS1[rows] = POS[rows], a[rows]);
				if var <> NULL and var >= 0 and type(var, integer) then
					RETURN(true);
				fi:
			fi:
		fi:
	od:

	RETURN(false);
end:

##############################################################################
# Plays Chomp starting at position POS
# and the human is first-to-move. 
#
# POS a solid position
##############################################################################

PlayChomp := proc(POS)
	local POS1, POS2, i, j, cell, randP, turn, len:
	
	for i from 2 to nops(POS) do
		if POS[i] > POS[i - 1] then
			printf(`You should have a non-increasing list of \n`):
			printf(`positive integers, e.g. PlayChomp([4,4,3]); \n`):
			ERROR(`Bad position`):
		fi;
	od:

	POS1 := POS:
	printf(`The starting position is \n`):
	DrawPos(POS1):
	
	while POS1 <> [1] do
		turn := 0;
		
		printf(`It is your turn \n`):

		printf(`Which Row?, type an integer between 1 and %d followed by ;<CR>. Type in q; to quit:\n`, nops(POS1)):

		j := readstat():
		if not type(j, integer) then
			if j = `q` then
				RETURN(`The game is terminated at user's request`);
			fi:
		else
			while not (1 <= j and j <= nops(POS1)) do 
				printf(`You must pick an integer between %d and %d\n`, 1, nops(POS1)):
				printf(`Try again!\n`):
				j := readstat():
			od:
		fi;
		
#		len := sum(POS1[ii], ii=1..j);
		len := POS1[j];
		printf(`Which Column? Type an integer between 1 and %d followed by <CR>. Type in q; to quit:\n`, len):

		i := readstat():
		if not type(i, integer) then
			if i = `q` then
				RETURN(`The game is terminated at user's request`);
			fi:
		else
			while not (1 <= i and i <= len) do 
				printf(`You must pick an integer between 1 and %d. Type in q; to quit:\n`, len):
				printf(`Try again!\n`):
				i := readstat():
			od:
		fi:
		
		cell := [i, j]:

		if cell = [1, 1] then
			printf(`Now there is nothing left, you lost!\n`):
			RETURN():
		fi:

		POS1 := Chop(POS1, cell):

		printf(`After your move, the new position is %a\n`, POS1):
		DrawPos(POS1):
		
		if POS1 = [1] then
			break;
		fi:
		
		turn := 1;
		
		POS2 := FindNextMoveWithNim(ConvertPosFromGameToProgFormat(POS1), 0, a):
		if POS2 = NULL then
			randP := rand(2..nops(POS1));
			i := randP();
			
			for j from i + 1 to nops(POS1) do
				if POS1[j] = POS1[i] then
					POS1 := [op(1..j-1, POS1), POS1[j] - 1, op(j+1..nops(POS1), POS1)];
				else
					break;
				fi:
			od:
			POS1 := [op(1..i-1, POS1), POS1[i] - 1, op(i+1..nops(POS1), POS1)];
			
			for i from 1 to nops(POS1) do
				if POS1[i] = 0 then
					POS1 := [op(1..i-1, POS1)];
					break;
				fi:
			od:
		else
			POS1 := ConvertPosFromProgToGameFormat(POS2);
		fi:

		printf(`I just moved. After my move, the position is %a\n`, POS1):
		DrawPos(POS1):
		
		if POS1 = [1] then
			break;
		fi:
	od:

	if POS1 = [1] then
		if turn = 1 then
			printf(`I won!, better luck next time.\n`):
		else
			printf(`You won!, Good job.\n`):
		fi:
		
		RETURN():
	fi:

end:

##############################################################################
# Given a position posit (a list of non-increasing
# positive integers), finds the result of perfoming the
# move of Chomping at the cell cell
#
# POS a solid position
# cell a two element list, the first is the column, and second row
##############################################################################

Chop := proc(POS, cell) 
	local i, j, POS1, j1: 
	
	i := cell[1]: 
	j := cell[2]:
	
	POS1 := [op(1..j-1, POS)]:
	for j1 from j to nops(POS) while POS[j1] >= i do
		if i > 1 then 
			POS1 := [op(POS1), i-1]: 
		fi:
	od:
	
	[op(POS1),op(j1..nops(POS),POS)]: 
end:

##############################################################################
# draws the position POS as a line of Xs
#
# POS a solid position
##############################################################################

DrawPos := proc(POS)
	local i, j, len:

	len := 0;
	for i from 1 to nops(POS) do
#		len := len + POS[i];
		len := POS[i];
		printf(`%s\n`, cat(seq(`X `, j = 1..len))):
	od:

	printf(`\n`):
end:

##############################################################################
# Convert a position POS from the format in the program to the format in the game
#
# POS a solid position
##############################################################################

ConvertPosFromProgToGameFormat := proc(POS)
	local i, j, ret;
	
	ret := [];
	
	for i from nops(POS) by -1 to 1 do
		ret := [op(ret), sum(POS[j], j = 1..i)]:
	od:
	
	ret:
end:

##############################################################################
# Convert a position POS from the format in the game to the format in the program
#
# POS a solid position
##############################################################################

ConvertPosFromGameToProgFormat := proc(POS)
	local i, j, ret;
	
	ret := [POS[nops(POS)]];
	
	for i from nops(POS) - 1 by -1 to 1 do
		ret := [op(ret), - POS[i + 1] + POS[i]]:
	od:
	
	ret:
end: