##############################################################################
# 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 P-position is previous player win, i.e. the nim value of the position equals to 0
#
# InitP is a list of P-positions that InitP[k][h] is a list
# of P-positions with k rows and the value of the first row equals to h
#
# this program returns all P-positions
##############################################################################

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

InitP := 'InitP':
Symb_Var := 'x': #'Symb_Var':
InitP := [[[[1]]], [[[1 + a[1], 1]]]]:

CHECK_VALUE := 1000000:
NON_FIXED_ROWS := 2:
MIN_OCCURANCE_NUMBER := 3:
PRINT_RESULT := false:

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 all chomp P-positions (privous player wins)`);
		printf(`%s\n`, `The main function is Main`);
		printf(`%s\n`, ``);
		printf(`%s\n`, `A chomp position is an ordered list`);
		printf(`%s\n`, `the number of operand is the same as 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`, `each of the other operands is the the number of excessive pieces in `);
		printf(`%s\n`, `the corresponding row. For example [1, 0, 1, 0, 2, 2] represents a `);
		printf(`%s\n`, `positions with 6 rows and has 1, 1, 2, 2, 4, 6 pieces in each row`);
		printf(`%s\n`, ``);
		printf(`%s\n`, `Help is available for the following functions:`);
		printf(`%s\n`, `Constants, Main, FindSymbolicWinningMove, FindWinningMove, 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 InitP or the empty list [] if it is the first time`);
			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 first row i`);
			printf(`%s\n`, `So the number of operands of the input/out is the same as the `);
			printf(`%s\n`, `number of maximum rows. In other words, T[i][j] 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 InitP = [[[[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 lager than 0`);
			printf(`%s\n`, ``);
			printf(`%s\n`, `Users can plug in a pre-calculated list T as the starting point.`);
		elif op(1, [args]) = FindSymbolicWinningMove then
			printf(`%s\n`, `The function takes two perameters: POS and a`);
			printf(`%s\n`, `where POS is the chomp position to find a winning move (P-position) for`);
			printf(`%s\n`, `a is a symbolic integer that the formulas will use`);
			printf(`%s\n`, `the return value is the P-position with possibly symbol a in it`);
		elif op(1, [args]) = FindWinningMove then
			printf(`%s\n`, `The function takes two perameters: POS and a`);
			printf(`%s\n`, `where POS is the chomp position to find a winning move (P-position) for`);
			printf(`%s\n`, `a is a symbolic integer that the formulas will use`);
			printf(`%s\n`, `unlike FindSymbolicWinningMove, the return value is an explicit P-position`);
			printf(`%s\n`, `(no symbols)`);
		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`, `InitP = `, InitP);
			printf(`%s\n`, ``);
			printf(`%s\n`, `InitP is the default P-position list. See Main function for detail`);
			printf(`%s\n`, ``);
			printf(`%s%d;\n`, `MIN_OCCURANCE_NUMBER = `, MIN_OCCURANCE_NUMBER);
			printf(`%s\n`, ``);
			printf(`%s\n`, `MIN_OCCURANCE_NUMBER is the minimal occurance of a pattern to be `);
			printf(`%s\n`, `considered as a formula. The default value of MIN_OCCURANCE_NUMBER `);
			printf(`%s\n`, `is 3, 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 is `);
			printf(`%s\n`, `enough to be considered as 24 + 2a when MIN_OCCURANCE_NUMBER = 3`);
			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 or 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:

##############################################################################
# Return all P-positions with max_rows rows and top (max_rows-2) rows have at most max_cols columns
#
# T[k][h] is a list of P-positions with k rows and the first rows of which are h when k > 2
# max_rows is the maximum number of rows to be calculated
# max_cols is the maximum number of cols in the top (max_rows - 2) rows
# a is a list of symbolic numbers
# returns a list of P-positions with the same format as T
##############################################################################

Main := proc(T, max_rows, max_cols, a)
	local i, j, len, LIST, max_deg, TT, retLIST;

	if nops(T) <= 2 then	# if T is "less than" InitP, use InitP instead
		TT := InitP;
	else
		TT := T;
	fi:
	
	if max_rows <= 2 then	# there is nothing to calculate
		RETURN(TT);
	else
		LIST := TT:
		len := nops(LIST):
		
		for i from 3 to max_rows do	# extend one row at a time
print(`CalculatePPosition: calculating less rows -- `, i, j):
			LIST := ExtendListOfPPositions(LIST, i, max_cols, a):
		od:
	fi:

print(`Formatting output. Please Wait...`);
	# reformat the result in sequencial order, e.g. [1,2,3] < [1,2,4]
	retLIST := [];
	for i from 1 to nops(LIST) do
		for j from 1 to nops(LIST[i]) do
			if j = 1 then
				retLIST := [op(1..i - 1, retLIST), [sort(LIST[i][j], SortChomp)]]:
			else
				retLIST := [op(1..i - 1, retLIST), 
					[op(1..j - 1, retLIST[i]), sort(LIST[i][j], SortChomp)]]:
			fi:
		od;
	od;
	
	# return the formatted list
	retLIST:
end:

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

PrintResult := proc(msg)
	if PRINT_RESULT or printlevel > 1 then
		print(msg);
	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:

#####################################################################################
# This function is to return the set of all the positions grown from POS 
# by adding new rows. 
# So if POS is P-position, then the result is the set of N-positions grown by 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):		
end:

#####################################################################################
# POS is a chomp position, grow POS to POS1 so that POS1[0] = POS[0] + c
# and POS can be deduced from POS1 by one move
# which means if POS is P-position, then the result is the set of 
# N-positions grown by POS with extra rows at the top
# NOTE: POS itself might be in the result
#
# 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:
#print(`GrowPositionFromTopByMaxValue: POS = `, POS, `c = `, c):

    if 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
				        RETURN( 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) ):
				    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
				ERROR(`calculation error`);
			fi:
        else	# c > POS[2] and POS[2] <> 0, cannot grow the position
            RETURN( {} ):
        fi:

    fi:

	# now POS[2] is symbolic
    if coeff(POS[2], a[2], 0) >= c then		# same as the instance above
        RETURN({[POS[1] + c, POS[2] - c, POS[3]]}):
    elif coeff(POS[2], a[2], 0) < c and coeff(POS[2], a[2], 0) > 0 then	# expand beyond the symbol
        RETURN({[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)), POS[3]]}):
    else # POS[2] = alpha * a[2], two cases: treat the symbol as a zero or not
        RETURN({[POS[1] + c, seq(POS[i], i = 2..nops(POS))]} 
			union GrowPositionFromTopByMaxValue([POS[1], 0, seq(POS[i], i = 3..nops(POS))], c, max_cols, a) ):
    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 P-position, then the result is the set of 
# N-positions grown by POS with the same number of rows and the same first row
# NOTE: POS itself might be in the result
#
# 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, list1, list2:
    
    list3 := {}: # the return value
    len := nops(POS):
#print(`GrowPositionFromBottomByMaxValue: POS = `, POS):
    
    for i from 2 to len - 1 do
        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
			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):
#print(`GrowPositionFromBottomByMaxValue: list1 = `, list1);
			
			# Put the value of the current row back into the list
			list1 := {seq( [ op(list2), 
				POS[i] + op(i, list1[jj]) 
					- iquo(coeff(POS[i] + op(i, list1[jj]), a[i], 1), 2) * a[i], 
					seq(list1[jj][j], j=i+1..len)], jj = 1..nops(list1)) }:

			list3 := list3 union list1:
			
			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 P-position
				# 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:

	# There is a possibility that the current position is in the list
    if type(POS[len], integer) and POS[len] = 0 then # do we really care????, it is going to happen any way
		list3 minus {POS}:
    else # the last line shall be expanded as symbolic, always
		(list3 union {[seq(POS[i], i=1..len - 1), 
		POS[len] + 1 + (1 - coeff(POS[len], a[len], 1)) * a[len]] }) 
		minus {POS}:
	fi:
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:

#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
		if rows <= 2 then
			RETURN( { [op(prefix), seq(a[nops(prefix) + i], i =1..rows)] } ):
		else
			RETURN( ExpandBlockByMaxValue(prefix, max_value, [a[nops(prefix) + rows - 1], 
				a[nops(prefix) + rows]], rows - 2, a) ):
		fi:
    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; #max_value * rows:
    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 = 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
			if nops(suffix) = 0 then
				ERROR(`should not be here`);
				# the last two expanded rows to infinity
		        RETURN( {seq([op(prefix), i, 
					a[nops(prefix) + 2] + sum_to_value - i], i = 0..sum_to_value - 1)}
			        union {[op(prefix), max(sum_to_value, 0) + a[nops(prefix)+1], a[nops(prefix) + 2]]} )
			else
				# 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)} ):
			fi:
        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:


##############################################################################
# returns two elements
# the first one is the generating funciton
# the second one is all the N-positions with the FixedTop
#
# T is the list of all P-positions
# 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, f, len, len1, POS, list1, list2, rows, firstrow;

	rows := nops(FixedTop) + 2;
	firstrow := FixedTop[1]:
	
#if FixedTop = [1, 1] then
#print(`BuildGeneratingFunctionForFixedTop: T = `, T, `  FixedTop = `, FixedTop, `  a = `, a);
#fi:

	f := 1:
	for i from rows - NON_FIXED_ROWS + 1 to rows do
		f := f / (1 - Symb_Var[i]):
	od:
#print(`BuildGeneratingFunctionForFixedTop: f = `, f);

	list1 := GetAllNPositionsWithFixedTop(T, FixedTop, max_cols, a):

	# calculate generating function by checking all N-positions
	list2 := list1[1];
	len := nops(list2):
	for i from 1 to len do:
		POS := list2[i]:
		f := f - ConvertPositionToGeneratingFunctionWithFixedTop(list2[i]);
	od:
#if FixedTop = [2, 6] then
#print(`BuildGeneratingFunctionForFixedTop: list2 = `, list2, `FixedTop = `, FixedTop, `f = `, f);
#fi:
	
	# calculate generating function by checking all P-positions
	list2 := list1[2];
	len := nops(list2):
	for i from 1 to len do:
		POS := list2[i]:
		f := f - ConvertPositionToGeneratingFunctionWithFixedTop(list2[i]);
	od:
#if FixedTop = [2, 6, 0] then
#print(`BuildGeneratingFunctionForFixedTop: list2 = `, list2, `T = `, T, `f = `, f, nops(list2));
#fi:
	
	f, list1[1]:
end:

##############################################################################
# Convert a position into generating function, 
# i.e. the monomial representation of the chomp position
#
# POS is a chomp position
##############################################################################

ConvertPositionToGeneratingFunctionWithFixedTop := proc(POS)
	local i, len, f;
	
#print(`ConvertPositionToGeneratingFunctionWithFixedTop: POS = `, POS, `f = `, f):

	f := 1:
	len := nops(POS);
	
	for i from len - NON_FIXED_ROWS + 1 to len do
		f := f * Symb_Var[i] ^ coeff(POS[i], a[i], 0):
		
		if coeff(POS[i], a[i], 1) <> 0 then
			f := f / (1 - Symb_Var[i] ^ coeff(POS[i], a[i], 1)):
		fi:
	od:
#print(`ConvertPositionToGeneratingFunctionWithFixedTop: POS = `, POS, `f = `, f):
	
	f:
end:

##############################################################################
# Increase the value of the first row of the computed P-positions of 
# the desired number of rows by 1
#
# T is a list of available P-positions
# rows is the number of rows we want to calculate
# max_cols is the maximum cols allowed for the top rows
# a is a list of symbolic variables
##############################################################################

ExtendListOfPPositions := proc(T, rows, max_cols, a)
	local i, j, jj, k, genfunc, newitems, len, newlen, bIntersect, 
		POS, POS1, NList, NList1, PTable, PTable1, listFixedTop, 
		FixedTop, lenFixedTop, AllNList;
	
	PTable := T:

	#calculate all possible top rows
	listFixedTop := [];
	newitems := ExpandBlockByMaxValue([], max_cols, 
		[max_cols + 1, max_cols + 1], rows - NON_FIXED_ROWS, a);
	len := nops(newitems);
	for i from 1 to len do
		if newitems[i][1] <> 0 and sum(newitems[i][ii], ii = 1..rows - NON_FIXED_ROWS) <= max_cols then
			listFixedTop := [ op(listFixedTop), [op(1..rows - NON_FIXED_ROWS, newitems[i])] ];
		fi:
	od:
	listFixedTop := sort(listFixedTop, SortChomp);
#print(`ExtendListOfPPositions: `, listFixedTop);
	
	for k from 1 to nops(listFixedTop) do
		FixedTop := listFixedTop[k];
		lenFixedTop := nops(FixedTop);
		
#print(`>>>>>>>>ExtendListOfPPositions<<<<<<<<: `, newitems, FixedTop);
		
		# For each top, recalculate the generating function again	
		AllNList := BuildGeneratingFunctionForFixedTop(PTable, FixedTop, max_cols, a);
		genfunc := AllNList[1]:
		AllNList := AllNList[2]:
		
#if FixedTop = [1, 1] then
#if nops(FixedTop) = 6 then
#print(`ExtendListOfPPositions: `, FixedTop, genfunc, mtaylor(genfunc, {x[8], x[7]}, 20));
#print(`ExtendListOfPPositions: `, FixedTop, PTable, genfunc);
#fi:

		# Find the next P-position with the top
		newitems := GetPositiveItemWithFixedTop(genfunc, FixedTop, Symb_Var);

#print(`ExtendListOfPPositions: `, PTable, listFixedTop, genfunc, newitems);

#		while nops(op(2, newitems)) > 0 do
		while nops(newitems) > 0 do
#print(`ExtendListOfPPositions: `, `Starting new column --- genfunc = `, genfunc, ` newitems = `, newitems):

#			newitems := op(2, newitems):
			newitems := newitems:
		
			PTable := AddNewPPositionIntoResult(PTable, newitems, rows, FixedTop[1], AllNList, a):
			newitems := PTable[1];
			newlen := nops(newitems);
			PTable := PTable[2];
		
#if FixedTop = [1, 1] then
#print(`ExtendListOfPPositions: AllNList = `, AllNList, ` genfunc = `, genfunc):
#fi:		

#print(`ExtendListOfPPositions: newitems = `, newitems, `  PTable[rows][FixedTop[1]] = `, PTable[rows][FixedTop[1]]):
		
			bIntersect := false;
			
			# Update the generating function to calculate the next P-position
			for i from 1 to newlen do			
				POS := newitems[i]:
				
				genfunc := genfunc - ConvertPositionToGeneratingFunctionWithFixedTop(POS):

				NList := convert(GrowPositionFromBottomByMaxValue(POS, max_cols, a), list):
				len := nops(NList);
				NList1 := [];
				for j from 1 to len do
					if [op(1..lenFixedTop, NList[j])] = FixedTop then
						NList1 := [ op(NList1), NList[j] ]:
					fi:
				od:
#print(`ExtendListOfPPositions: NList = `, NList, `  NList1 = `, NList1, `  FixedTop = `, FixedTop, len, [op(1..nops(FixedTop), NList[1])]):
				NList := [op(NList1), POS];	# add in POS to check all previous P-positions
				
				len := nops(PTable[rows][FixedTop[1]]);

				for j from 1 to len do
					POS1 := PTable[rows][FixedTop[1]][j];
					if POS1 <> POS 
							and not ElementIsASolution(POS, POS1, a)
							and ElementsIntersect(POS1, NList, a) then
PrintResult(cat(`!!!!! ExtendListOfPPositions: element REMOVED `, PTable[rows][FixedTop[1]][j], `by `, POS)):

						# remove the element because the pattern we guessed conflict with 
						# the newly found position
						# ????but the pattern will generate faulty P-positions
						# are they all going to be corrected eventually or not?
						# how many more positions shall we remove????
						if type(PTable[rows][FixedTop[1]][j][rows -1], integer) then
							# We are removing a solid position instead of a pattern
							# something must be wrong
							ERROR(`Removing a solid position. Something is wrong. Please check the 3-row base is large enough`);
						fi:
						
						# remove the pattern and add in MIN_OCCURANCE_NUMBER items in the pattern
						# to speed up the program (by a little) 
						PTable := [ op( 1..rows - 1, PTable ),
							[ op( 1..FixedTop[1] - 1, PTable[rows] ),
							sort( [op( 1..j - 1, PTable[rows][FixedTop[1]] ),
							seq( subs(a[rows - 1] = jj, PTable[rows][FixedTop[1]][j]), jj = 0..MIN_OCCURANCE_NUMBER - 1 ),
							op( j + 1..nops(PTable[rows][FixedTop[1]]), PTable[rows][FixedTop[1]] ) ], SortChomp),
							op( FixedTop[1] + 1..nops(PTable[rows]), PTable[rows] ) ],
							op( rows + 1..nops(PTable), PTable ) ];

						bIntersect := true;
					fi:
				od:
				
				if bIntersect then
					# Since something is removed from the P-list, we need to recalculate the
					# generating function
					AllNList := BuildGeneratingFunctionForFixedTop(PTable, FixedTop, max_cols, a);
					genfunc := AllNList[1]:
					AllNList := AllNList[2]:

					break;
				else
					# Since nothing wrong happened, simply add the N-positions 
					# into the gernerating funciton
					NList := NList1;
					AllNList := [ op(AllNList), op(NList) ]:
				
					len := nops(NList);
					for j from 1 to len do
						genfunc := genfunc - ConvertPositionToGeneratingFunctionWithFixedTop(NList[j]);
					od:
				fi:
			od:
			
			# Now find the next P-position
			newitems := GetPositiveItemWithFixedTop(genfunc, FixedTop, Symb_Var);
		od:
	od:
	
	PTable:
end:

##############################################################################
# Get the "first" weighted item with positive coeff in the numerator of 
# a quotation P(x1,...,xn)/(1-x1)...(1-xn)
# if an item is found, we have to make sure there is no items with less degree
# with "bigger" negative coeff due to the denominator, e.g. - x[2] + x[2]x[3] is
# considered all negative (check the taylor expansion!)
# return a list of items with lowest "weighted" degree, [] if nothing found
#
# genfunc is the function we are checking against
# FixedTop is the top rows of the chomp positions that genfunc was generated from
# a is a list of symbolic variables
##############################################################################

GetPositiveItemWithFixedTop := proc(genfunc, FixedTop, a)
	local i, j, len, list1, returnlist, size, numerf, denomf, tempf,
		itemdeg, itemcoeff, itemtotaldeg, maxdeg, returndeg, f;

	size := nops(FixedTop) + NON_FIXED_ROWS:
	f := genfunc;
	numerf := numer(f);
	denomf := denom(f);
#	for i from size - NON_FIXED_ROWS + 1 to size do
#		# calculate genfunc * (1 - a[i])
#		# seperate the denominator and numerator so that the polynomials are less complicated
#		tempf := denomf / (1 - a[i]):
#		if denom(tempf) = 1 then
#			denomf := tempf;
#		else
#			numerf := numerf * (1 - a[i]):
#		fi:	
#	od:
	numerf := expand(numerf);
	denomf := expand(denomf);
#print(`GetPositiveItemWithFixedTop: 2 `, genfunc, f, denom(f), numer(f), size, a):
	
	#???? how big shall be the maximum degree to be safe????
	maxdeg := 2;
	for i from 1 to size do
		maxdeg := maxdeg + degree(numerf, a[i]) + degree(denomf, a[i]):
	od;

#print(`GetPositiveItemWithFixedTop: `, f, {seq(a[i], i=size - NON_FIXED_ROWS + 1..size)}, maxdeg):
	f := mtaylor(genfunc, {seq(a[i], i=size - NON_FIXED_ROWS + 1..size)}, maxdeg):

#if FixedTop = [1, 0] then
#print(`GetPositiveItemWithFixedTop: FixedTop = `, FixedTop, `f = `, f, `size = `, size, `a = `, a):
#fi:
	
	list1 := GetItemsFromPolynomial(f):
	len := nops(list1):
	returnlist := [];
	returndeg := CHECK_VALUE;
	
	for i from 1 to len do
		itemcoeff := list1[i]:
		itemtotaldeg := 0:

		for j from size - NON_FIXED_ROWS + 1 to size do
			itemdeg[j] := degree(itemcoeff, a[j]);
			itemcoeff := coeff(itemcoeff, a[j], itemdeg[j]);
			itemtotaldeg := itemtotaldeg + itemdeg[j] * (size - j + 1):	#degree is "weighted"
		od:
#print(`GetPositiveItemWithFixedTop: itemcoeff = `, itemcoeff, `returnlist = `, returnlist):	
		
		if itemcoeff > 0 then
			if itemtotaldeg < returndeg then
				returndeg := itemtotaldeg:
				returnlist := [[op(FixedTop), seq(itemdeg[j], j=size - NON_FIXED_ROWS + 1..size)]];
			elif itemtotaldeg = returndeg then
				returnlist := [op(returnlist), 
					[op(FixedTop), seq(itemdeg[j], j=size - NON_FIXED_ROWS + 1..size)]];
			fi:
		fi:
	od;	
	
#print(`GetPositiveItemWithFixedTop: ii = `, ii, `returnlist = `, returnlist):	
#	[returndeg, returnlist]:
	returnlist:
end:

##############################################################################
# Get all the N-positions from P-positions list T
# all N-positions as the first return element, 
# all P-positions of the same row/column as the second return element
#
# T is a list of all available P-positions
# FixedTop is the value of the top rows
# max_cols is the maximum cols allowed for the top rows
# a is a list of symbolic variables
##############################################################################

GetAllNPositionsWithFixedTop := proc(T, FixedTop, max_cols, a)
	local i, j, k, l, len, len1, POS, listT, LIST, lenLIST, lenlistT, 
		returnset, returnset2, set1, rows, firstrow;
	
	rows := nops(FixedTop) + 2;
	firstrow := FixedTop[1];
	
#print(`GetAllNPositionsWithFixedTop: T = `, T, `  FixedTop = `, FixedTop, `  a = `, a);
	returnset := {};
	returnset2 := {};
	
#print(`GetAllNPositionsWithFixedTop: f = `, f);

	if nops(T) >= rows then
		len := nops(T[rows]):
		if len >= firstrow then
			len := firstrow - 1:
		fi:
#print(`GetAllNPositionsWithFixedTop: len = `, len, `firstrow = `, firstrow):
		
		# Get N-positions from all P positions with the same number of rows
		
		# check N-positions 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]:
				set1 := GrowPositionFromTopByMaxValue(POS, firstrow - i, max_cols, a):
				len1 := nops(set1);
				for k from 1 to len1 do
					if [op(1..rows - NON_FIXED_ROWS, set1[k])] = FixedTop then
						returnset := returnset union { set1[k] }:
					fi:
				od:
			od:
		od:
		
		# check N-positions with the first row the same as this one
		if nops(T[rows]) >= firstrow then
			LIST := T[rows][firstrow];
			lenLIST := nops(LIST);
			
			for j from 1 to lenLIST do
				if [op(1..rows - NON_FIXED_ROWS, LIST[j])] = FixedTop then
					returnset2 := returnset2 union { LIST[j] } :
				fi:
			od:
			
			for j from 1 to lenLIST do
				POS := LIST[j];
				set1 := GrowPositionFromBottomByMaxValue(POS, max_cols, a):
				len1 := nops(set1);
				for k from 1 to len1 do
					if [op(1..rows - NON_FIXED_ROWS, set1[k])] = FixedTop then
						returnset := returnset union { set1[k] }:
					fi:
				od:
			od:
		fi:
	fi:
#print(`GetAllNPositionsWithFixedTop: firstrow = `, firstrow, ` returnset = `, returnset);

	# Get N-positions from all P positions with less rows

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

		for j from 1 to len do
			len1 := nops(listT[j]);
			for k from 1 to len1 do
				set1 := ExpandPositionWithNewRow(listT[j][k], rows, firstrow, max_cols, a):

				for l from 1 to nops(set1) do
					if	[op(1..rows - NON_FIXED_ROWS, set1[l])] = FixedTop then
						returnset := returnset union { set1[l] }:
					fi:
				od:
			od:
#print(`GetAllNPositionsWithFixedTop: listT = `, listT, `j = `, j, `LIST = `, LIST):

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

	[returnset, returnset2];	
end:

##############################################################################
# Check if a position POS 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
#
# POS is a position
# LIST is a list of positions
# 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):
	POS1 := subs(a = aa, POS):
	
	bIntersect := false:
#print(`ElementsIntersect: POS1 = `, POS1, `LIST = `, LIST):
	
	for i from 1 to len do
		POS2 := LIST[i];
		
		bIntersect := true:
#print(`ElementsIntersect: POS1 = `, POS1, `POS2 = `, POS2, `res = `, res):
			
		for j from 2 to rows do
			var := POS1[j] - POS2[j];

			if type(var, numeric) then
				if var <> 0 then
					bIntersect := false:
					break;
				else
					next:
				fi:
			else
				coefaa := coeff(var, aa[j], 1);
				coefa := coeff(var, a[j], 1):
				const := coeff(coeff(var, aa[j], 0), a[j], 0);
				coef := gcd(coefa, coefaa);
									
#print(`ElementsIntersect: POS1 = `, POS1, `POS2 = `, POS2, `const = `, const, `coef = `, coef, `j = `, j);
				# we need an positive integral solution
				if type(const / coef, integer) then
					if (coefaa = 0 and const / coefa > 0) or (coefa = 0 and const / coefaa > 0) then
							
						# solvable, but only with negative solution
						bIntersect := false:
						break;
					fi:
				else	
					# solvable, but not integral solution
					bIntersect := false:
					break;
				fi:
			fi:
		od:
#print(`ElementsIntersect: POS1 = `, POS1, `POS2 = `, POS2, `res = `, res, `bIntersect = `, bIntersect):
		
		if bIntersect then
#print(`ElementsIntersect: POS1 = `, POS1, `POS2 = `, POS2, `bIntersect = `, bIntersect):
			break;
		fi:
	od:
	
	bIntersect:
end:

##############################################################################
# Check if POS2 is a special case of POS1
# if both POS1 and POS2 is symbolic, return false
#
# POS1 is a position
# POS2 is a position
# a is a list of symbolic variables
##############################################################################

ElementIsASolution := proc(POS1, POS2, a)
	local i, bSymbolic, rows;
	
	bSymbolic := false;
	rows := nops(POS1);
	
	for i from rows - NON_FIXED_ROWS + 1 to rows do
		if coeff(POS1[i], a[i], 1) <> 0 then
			bSymbolic := true;
			break;
		fi;
	od:
	
	if bSymbolic then
		for i from rows - NON_FIXED_ROWS + 1 to rows do
			if coeff(POS2[i], a[i], 1) <> 0 then
				RETURN(false);
			fi;
		od:
	fi:
	
	RETURN( ElementsIntersect(POS1, [POS2], a) );
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 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:

##############################################################################
# Add new positions into a result P-position 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 P-positions
# the second is the list of the new all P-positions
#
# PTable is the list of P-positions
# value is a list of P-positions 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 N-positions
# a is a list of symbolic variables
##############################################################################

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

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

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

##############################################################################
# Given a solid P-position, try to find out if there is a pattern forming in 
# the P-positions PList. A pattern is when certain criteria happens at least
# MIN_OCCURANCE_NUMBER times, and the pattern does not conflict with any
# of the N-positions 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 solid P-position
# PList is the list of all P-positions
# NList is the list of all N-positions
# a is a list of symbolic variables
##############################################################################

FindPattern := proc(POSAdd, PList, NList, a)
	local len, gap, firstoccuranceat, occurance, firstindex, i, POS, rows, lastgap, patternlist;
	
	len := nops(PList);
	rows := nops(POSAdd);
	
	lastgap := 0;
	patternlist := [];

	while true do	# get all the possible patterns available
		firstindex := POSAdd[rows - 1];
		gap := 0:
		firstoccuranceat := len;
		occurance := 1:

		# find the first index of the entry after which all the entries match the proposed formula
		for i from len 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 P-position that might generate a symbolic sequence
				# firstindex is the value of the (rows - 1)th row of the first occurence
				# firstoccuranceat is the index of the first occurence in the list PList
				# occurance 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 
						POS[rows] = POSAdd[rows] and
						POSAdd[rows - 1] - POS[rows - 1] > lastgap then
					
					firstindex := POS[rows - 1]:
					firstoccuranceat := i;
					occurance := occurance + 1:
					gap := POSAdd[rows - 1] - POS[rows - 1]:
#print(` FindPattern creating gap: `, ` POS = `, POS, ` firstindex = `, firstindex, ` firstoccuranceat = `, firstoccuranceat, ` occurance = `, occurance, ` gap = `, gap);
				fi:
			else
				# find the first possible P-position that might generate a pattern
				if op(1..rows - 2, POS) = op(1..rows - 2, POSAdd) and 
						firstindex - POS[rows - 1] = gap then
					if POS[rows] = POSAdd[rows] then
						firstindex := POS[rows - 1]:
						occurance := occurance + 1:
								
#print(` FindPattern increasing occurance: `, ` POS = `, POS, ` firstindex = `, firstindex, ` firstoccuranceat = `, firstoccuranceat, ` occurance = `, occurance, ` gap = `, gap);
						if occurance >= MIN_OCCURANCE_NUMBER then
							break;
						fi:
					else
#print(` FindPattern elminating gap: `, ` POS = `, POS, ` firstindex = `, firstindex, ` firstoccuranceat = `, firstoccuranceat, ` occurance = `, occurance, ` gap = `, gap);
						gap := 0;
						occurance := 1;
						i := firstoccuranceat:
					fi:
				fi:
			fi:
			
			# the pattern is invalid since we are unable to find the next desired position
			# loop back to the previous spot
			if i = 1 and gap <> 0 and firstoccuranceat <> 1 then
				gap := 0;
				occurance := 1;
				i := firstoccuranceat;
			fi:
		od:
		
		if gap <> 0 and  occurance >= MIN_OCCURANCE_NUMBER then
			lastgap := gap;
			patternlist := [op(patternlist), [op(1..rows - 2, POSAdd), firstindex + gap * a[rows - 1], POSAdd[rows]]];
		else
			break;
		fi:
	od:
#print(` FindPattern found patterns: `, patternlist);
				
	len := nops(patternlist);
	if len > 0 then

		for i from 1 to len do
			POS := patternlist[i];
			
#print(` find enough times: `, ` POS = `, POS, ` firstindex = `, firstindex, ` firstoccuranceat = `, firstoccuranceat, ` occurance = `, occurance, ` gap = `, gap);
			if not ElementsIntersect(POS, NList, a) then
PrintResult(cat(` Found SYMBOLIC position: `, ` POS = `, POS));
				RETURN( POS, sort( convert(convert(PList, set) minus 
					{ seq( subs(a[rows- 1] = i, POS), 
					i = 0..MIN_OCCURANCE_NUMBER - 1) } union { POS }, list), SortChomp) );
			fi:
		od:

PrintResult(cat(` Find position INTERSECTS: `, ` POSAdd = `, POSAdd)):
		RETURN(	POSAdd, sort( [op(PList), POSAdd], SortChomp ) );
	else
PrintResult(cat(` Find SINGLE position: `, ` POSAdd = `, POSAdd)):
		POSAdd, sort( [op(PList), POSAdd], SortChomp);
	fi:
end:

##############################################################################
# Given a posotion, find a solid winning move
#
# POS is a chomp position
# a is a list of symbolic variables
##############################################################################

FindWinningMove := proc(POS, a)
	local POS1, POS2, i, rows, rowsize, rowsize1, var;
	
	rows := nops(POS);
	POS1 := FindSymbolicWinningMove(POS, a);
	
	if POS1 = NULL then
		PrintResult(cat(`This is a P-position: `, POS));
		RETURN();
	elif POS1 = [1] then
		RETURN([1]);
	elif nops(POS1) < rows then
		[sum(POS[ii], ii = 1 .. rows - nops(POS1)), op(rows - nops(POS1) + 2 .. rows, POS)];
	else		
		POS2 := [op(1..rows - NON_FIXED_ROWS, POS1)];

#Only works for two non-fixed rows
		rowsize := sum(POS[ii], ii=1..rows):
		rowsize1 := sum(POS1[ii], ii=1..rows):
		
		var := solve(rowsize - rowsize1 = 0, {seq(a[i], i=(rows - NON_FIXED_ROWS + 1)..rows)});
		for i from rows - NON_FIXED_ROWS + 1 to rows do
			if subs(a[i] = subs(var, a[i]), POS1[i]) > POS[i] then
				POS2 := [op(POS2), subs(a[i] = 0, POS1[i])];
			else
				POS2 := [op(POS2), subs(a[i] = subs(var, a[i]), POS1[i])];
			fi:
		od:
			
		POS2:
	fi:
end:

##############################################################################
# Given a position, find a winning move. The result can be symbolic
#
# POS is a chomp position
# a is a list of symbolic variables
##############################################################################

FindSymbolicWinningMove := proc(POS, a)
	local i, j, k1, k2, k3, len1, len2, len3, rows, NList, NList1, FixedTop, 
		eq, res, var, POS2, max_cols;
	
	rows := nops(POS):
	FixedTop := [op(1..rows - NON_FIXED_ROWS, POS)];
	
	if rows > nops(InitP) then
		ERROR(`Cannot find winning move: POS beyond the limit of the result`);
	elif rows = 1 then
		if POS[1] > 1 then
			RETURN( [1] );
		fi;	
	elif rows = 2 then
		if POS[2] > 1 then
			RETURN( [POS[1], 1] );
		elif POS[2] = 0 then
			if POS[1] = 1 then
				RETURN( [1] );
			else
				RETURN( [POS[1] - 1, 1] );
			fi:
		fi:
	else
		max_cols := sum(POS[ii], ii = 1..rows - NON_FIXED_ROWS);
		
		for k1 from 1 to rows - 1 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];
					NList := ExpandPositionWithNewRow(POS2, rows, POS[1], max_cols, a);
#print(`FindWinningMove 1: POS2 = `, POS2, ` NList = `, NList);

					if ElementsIntersect(POS, NList, a) then
						RETURN( POS2 );
					fi:
				od;
			od:
		od:
		
		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];
				NList := GrowPositionFromTopByMaxValue(POS2, POS[1] - POS2[1], max_cols, a);
#print(`FindWinningMove 2: POS2 = `, POS2, ` NList = `, NList);
			
				if ElementsIntersect(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];
				NList := GrowPositionFromBottomByMaxValue(POS2, max_cols, a);
#print(`FindWinningMove 3: POS2 = `, POS2, ` NList = `, NList);
			
				if ElementsIntersect(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((`FindWinningMove: Nothing found for `, POS));
	RETURN():
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:

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

	printf(`\n`):
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:
	
	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>\n`, nops(POS1)):

		j := readstat():
		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:
		
		printf(`Which Column? Type an integer between 1 and %d followed by <CR>\n`, POS1[j]):

		i := readstat():
		while not (1 <= i and i <= POS1[j]) do 
			printf(`You must pick an integer between 1 and %d\n`, POS1[j]):
			printf(`Try again!\n`):
			i := readstat():
		od:

		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 := FindWinningMove(ConvertPosFromGameToProgFormat(POS1), 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:

##############################################################################
# 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: