#read `result_C.txt`:
readlib(mtaylor):

# each position (a, b, c) has a beans in the first pile, a+b beans second, and a+b+c third
#
#

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

Help := proc()
	if args = NULL then
		print(`This package verifies the conjecture for N-heap Wythoff's game`);
		print(``);
		print(`N-heap Wythoff's game is a two-player impartial game`);
		print(`consisting of n piles of beans`);
		print(`Each player take turns to take any number of beans from any one pile,`);
		print(`or he can take (a1, ..., an) beans from all piles correspondingly, `);
		print(`providing the nim sum of the numbers is 0.`);
		print(``);
		print(`There are two conjectures involving the game. Please refer to`);
		print(`the associated paper "On Fraekel's N-heap Wythoff's Conjecture" for details.`);
		print(`This package proves that the conjectures are correct for the games with three-heap`);
		print(``);
		print(`The following is a list of functions provided in the package:`);
		print(`Main, Help`);
		print(``);
		print(`Type Help(functionname) for detailed help for each function.`);
		print(`Please check at www.math.temple.edu/~xysun for the details of the algorithm,`);
		print(`and report any questions or concerns to `);
		print(`xysun@math.temple.edu or zeilberg@math.rutgers.edu`);
	elif args[1] = Help then
		Help();
	elif args[1] = Main then
		print(`Main function takes one parameter which is the maximum size for the first heap,`);
		print(``);
		print(`It returns the following data in order: P, T, S, $\\alpha$, [N1, N2, N3]`);
		print(`the usage of the variables are same as described in the paper.`);
		print(``);
		print(`For example, Main(5) verifies the conjectures for the three-heap games`);
		print(`with the first heap having up to 5 beans.`);
	fi;
end:

# Given a three-piled P-position (convertfrom, x, y), find the gernerating functions for 
# all of its induced three-piled N-positions with convertto stones in the first pile.
ConvertOnePrevPPos := proc(convertfrom, convertto, x, y, x1, x2)
	local i, disp, conv, xcoord, ycoord, f, mp;

	disp := NULL;
	conv := convertto - convertfrom;
	f := 0;
	
	# the nim sum is of period mp
	mp := 2^(floor(simplify(log[2](conv))) + 1);
	
	# all new positions obtained by adding (a, b, c) whose nim sum is 0 are N
	for i from 0 to mp - 1 do
		xcoord := x - conv + i;
		ycoord := y + NimSum(i, conv) - i;

		# if xcoord or ycoord < 0, increase them
		if ycoord < 0 then
			xcoord := ycoord + xcoord:
			ycoord := - ycoord:
		fi;
		
		if xcoord < 0 then
			xcoord := xcoord + mp * ceil(-xcoord / mp):
		fi;
		
		f := f + x1^xcoord * x2^ycoord / (1 - x1^mp);
	od;
	
	# simply adding pieces to the first pile creates N
	if x >= conv then
		f := f + x1^(x - conv) * x2^y;
	fi;
	
	# add same number of pieces to the first and third piles
	if x >= conv then
		f := f + x1^(x - conv) * x2^(y + conv);
	fi;
	
	# add same number of pieces in the first and second piles
	if conv > y then
		if y + x - conv >= 0 then
			f := f + x1^(y + x - conv) * x2^(- y + conv);
		fi;
	elif x >= 0 then
		f := f + x1^x * x2^(y - conv);
	fi;

	f;
end:

# List all the P-positions whose first pile is less then currindex 
# and mark all their induced three-piled N-positions
# new_x1 and new_x2 is the possible new position
# and deg is the maximum value of x1 + x2 to be calculated.
# x1 and x2 are symbolic variables for the second and third number is our notation.
# return the generating function for all the relevant N-positions
ConvertAllPrevPPos := proc(PPos, currindex, new_x1, new_x2, deg, x1, x2)
	local i, j, NIndex, x, y, f;
	
	f := 0;

	for i from 0 to currindex - 1 do
		for j from 1 to nops(PPos[i]) do
			x := PPos[i][j][1];
			y := PPos[i][j][2];
			
			# Exclude prev P-pos that are not going to affect current calculation
			if (y < new_x2 - (currindex - i)) or (y > currindex - i + deg - new_x1) then
				next;
			fi;
#			if (i + x < currindex + new_x1 and y < new_x2 - (currindex - i))
#					or (i + x > currindex + max_x1) then
#				next;
#			fi;
			
			# convert current positions
			f := f + ConvertOnePrevPPos(i, currindex, x, y, x1, x2);
				 	
			# if the second pile is less than currindex, move it to the first and convert
		 	if i + x < currindex then
				f := f + ConvertOnePrevPPos(i + x, currindex, - x, x + y, x1, x2);
		 	fi;
		 	
			# if the third pile is less than currindex, move it to the first and convert
		 	if i + x + y < currindex then
				f := f + ConvertOnePrevPPos(i + x + y, currindex, - x - y, x, x1, x2);
		 	fi;		 	
		od:
	od;

	f;
end:

# Given a P-position, find all of its induced N-positions whose first pile remains the same
# mex_X, mex_Y are the lower bound of the degrees of all relevant monomials
# deg is the upper bound of the degrees of all relevant monomials
# returns the generating function of all the relevant parts
ConvertOneCurrentPPos := proc(x, y, mex_X, mex_Y, deg, x1, x2)
	local i, f;

	# if we have a P-position
	# we cannot have another P-position with the same number of pieces 
	# in the second or the third pile,
	# nor the same differences between the piles.
	
#	f := x1^x * x2^y / (1 - x2) 
#		+ x1^x * sum(x2^ii * x1^(y -ii), ii=0..y) + x1^(x+y) / (1 - x2)
#		+ x1^x * x2^y / (1 - x1);

	f := 0;
	
	if x >= mex_X then
		f := f + x1^x * x2^y / (1 - x2);
	fi;
	
	if y >= mex_Y then
		f := f + x1^x * x2^y / (1 - x1);
	fi;
	
	if x + y >= mex_X and x + y <= deg then
		f := f + x1^(x+y) / (1 - x2);
	fi;
	
	for i from 0 to y do
		if x + y - i >= mex_X and i >= mex_Y and x + y - i <= deg then
			f := f + x2^i * x1^(x + y -i);
		fi;
	od;
	
	f:
end:

# Find the generating functions of all special cases for three-piled (P and N) positions 
# with index pieces in the first pile
# index is the number of pieces in the first pile
# PPos is the list of all P-positions
# returns the generating functions and the numbers X must skip
ConvertSpecialPrevPPos := proc(PPos, index, x1, x2)
	local i, j, disp, coordx, f, skip;

	f := 0;
	skip := [];
	
	# go through previous P-positions to look for special cases 
	# not included in function ConvertAllPrevPPos
	for i from 0 to index-1 do
		for j from 1 to nops(PPos[i]) do
			# if the second pile has index pieces, 
			# we can add any number of pieces to the first pile to create N-positions
			# by moving the second pile to the first, the new second pile will grow from 0 to 
			# the size of the third, and then it becomes the third pile, which goes to infinity.
			# OR, we can add pieces to the first and the third piles at the same time
			if i + PPos[i][j][1] = index then
				f := f + sum(x2^ii * x1^(PPos[i][j][2] -ii), ii=0..PPos[i][j][2])
						+ x1^PPos[i][j][2] / (1 - x2)
						+ x2^(PPos[i][j][1] + PPos[i][j][2]) / (1 - x1);
				skip := [op(skip), PPos[i][j][2]];
			fi;

			# if the third pile has index pieces,
			# we can add pieces to the first and the second at the same time
			if i + PPos[i][j][1] + PPos[i][j][2] = index then
				f := f + x2^PPos[i][j][1] / (1 - x1);
			fi;
		od;
	od;
	
	f, skip;
end:

# Find P-positions whose first piles have up to num pieces
# returns in order: P-positions, numbers X must skip, numbers Y must skip, 
# differences between X and Y * golden_section, start of patterns
Main := proc(num)
	local i, j, l, SkipX, SkipY, f1, g_skipY, g, g1, deg, k1, k2, x, y, 
		mex_X, mex_Y, prev_y, indx_y, prev_diff, indx_diff, tempPPos, 
		golden, diffY, diff_golden, X, temp, res, checkholes, holes_filled,
		extend_index, f_extend, f_skipY, 
		formalargn, realargn, i_period, pattern_start,
		PPos1, SkipX1, SkipY1, DiffY1, PatternIndx1,
		MAX_EXT_SIZE, MAX_DEG_INCREMENT, MIN_CONVERT_EXT, EXPANSION_SIZE, MIN_DIFFY_DIFF;

	# constants	
	MAX_EXT_SIZE := 3;
	MIN_DIFFY_DIFF := 4:

	indx_y := -10;	# the starting index such that the third #s are incremented by 1
	prev_y := -10;	# the previous third #
	golden := (1 + sqrt(5)) / 2;	# the golden section

	# check input for previously calculated results
	formalargn := nops([op(1, op(procname))]);
	realargn := nops([args]);
	if realargn = formalargn + 5 then
		PPos1 := args[formalargn + 1];
		SkipX1 := args[formalargn + 2];
		SkipY1 := args[formalargn + 3];
		DiffY1 := args[formalargn + 4];
		PatternIndx1 := args[formalargn + 5];
	fi;

	# start the loop
	for i from 0 to num do
		# if the previous P-positions are given, skip i
		if type(PPos1[i], list) then 
			if type(PatternIndx1[i], list) and 
					nops(PPos1[i]) > PatternIndx1[i][2] then
				next;
			fi;
		else		
			PPos1[i] := [];
		fi;

		# find the optimal increment of deg
		MAX_DEG_INCREMENT := 30;
		if i > 0 then
			MAX_DEG_INCREMENT := max(MAX_DEG_INCREMENT, 
				floor((PPos1[i-1][PatternIndx1[i-1][2]][1] + 
				PPos1[i-1][PatternIndx1[i-1][2]][2]) / 8));
		fi;	
#print(i, `MAX_DEG_INCREMENT is `, MAX_DEG_INCREMENT, `time = `, time()):
		
		# (nim addition of i and a) minus a will be periodic of period i_period
		if i > 0 then
			i_period := 2^(floor(simplify(log[2](i))) + 1);
		else
			i_period := 2;
		fi;

		MIN_CONVERT_EXT := MAX_DEG_INCREMENT + i_period:
		EXPANSION_SIZE := MAX_DEG_INCREMENT;
		
		checkholes := false;	# indicates if certain numvers of y has been skipped
		pattern_start := 0;		# indicates the next position to test the pattern
		extend_index := -1;		# indicates the first index of the current run of 
								# finding P-positions without using the generating function

		# init variables
		if nops(PPos1[i]) > 0 then
			X := [];
			diffY := [];
			diff_golden := [];
			for l from 1 to nops(PPos1[i]) do
				x := PPos1[i][l][1];
				y := PPos1[i][l][2];
				X := [op(X), x, x + y];
				diff_golden := [op(diff_golden), i + x - floor(y * golden)];
				diffY := [op(diffY), y];
			od;
			
			if type(SkipX1[i], list) then
				X := [op(X), op(SkipX1[i])]:
				SkipX := SkipX1[i];
			else
				SkipX := [];
			fi;

			if type(SkipY1[i], list) then
				SkipY := SkipY1[i];
				diffY := [op(diffY), op(SkipY)];
			else
				SkipY := [];
			fi;			
		else
			res := ConvertSpecialPrevPPos(PPos1, i, x1, x2);
			X := res[2];		# all the numbers used by x and y in P-positions [i, x, y]
			diffY := [];		# all the numbers used by y
			diff_golden := [];	# the differences of i + x - floor(y * golden )
			SkipX := res[2];	# the numbers cannot occur in  x and x + y
			SkipY := [];		# the numbers cannot occur in y
		fi;

		g := 0;					# the generating function mtaylored to degree deg
		f_skipY := 0;			# the generating function to find SkipY mtaylored to degree deg
		deg := 0;				# max degree of mtayloring functions
		x := -1;				# [i, x, y] will be a P-position having i, i+x, i+x+y pieces
		y := -1;
		
		# look for the next P-position until a pattern is found, 
		# i.e.the conjectures are proved
		while true do
			mex_X := Mex(X);		# the next x value
			mex_Y := Mex(diffY);	# the minimum of the next y value
			holes_filled := false;	# indicated SkipY has been updated or not
			
			# find the next positions w\o using the generating function
			res := FindNextPPos(PPos1, i, mex_X, diffY, mex_Y, DiffY1, 
				PatternIndx1, MIN_DIFFY_DIFF);	
			
			if res = NULL then	# if we cannot find one
				if extend_index > 0 then	# update the generating function if necessary
					for l from extend_index to nops(PPos1[i]) do
						f1 := ConvertOneCurrentPPos(PPos1[i][l][1], PPos1[i][l][2], 
							mex_X, mex_Y, deg, x1, x2);
						g := g - mtaylor(f1, [x1, x2], deg);
					od;
					g := g - mtaylor(f_extend, [x1, x2], deg);
				else
					if x >= 0 then
						f1 := ConvertOneCurrentPPos(x, y, mex_X, mex_Y, deg, x1, x2);
						g := g - mtaylor(f1, [x1, x2], deg);
					fi;
				fi;
				
				x := -1;
				k1 := mex_X;			
				k2 := mex_Y;
				extend_index := -1;
				
				while x < 0 do
					g1 := coeff(g, x1, k1);

					for k2 from k2 to deg - k1 - 1 do
						if coeff(g1, x2, k2) > 0 then
							x := k1;
							y := k2;
							break;
						fi;
					od;

					# update the generating function by updating SkipY
					if x < 0 then
						# all holes in diffY must be filled
						if checkholes then
							res := TestSkippedNum(i, mex_X, max(op(diffY)), 
								diffY, f_skipY, deg, i_period, x1, x2);			
							checkholes := not res[1];
							SkipY := [op(SkipY), op(res[2])];
							diffY := [op(diffY), op(res[2])];
if res[2] <> [] then
print(i, `skipped Y: `, SkipY);
fi;

							if nops(res[2]) > 0 then
								holes_filled := true;
								break;
							fi;
						fi;
					fi;

					# update the generating function by increasing deg
					if x < 0 then
						deg := max(deg + MAX_DEG_INCREMENT, mex_X + max(op(diffY)) 
							+ MAX_DEG_INCREMENT) + i_period + 2;
						res := RecalculateGenFunc(PPos1, i, mex_X, mex_Y, diffY, deg, 
							MIN_CONVERT_EXT, checkholes, x1, x2);
					
						g := res[1];
						f_skipY := res[2];
					fi;
				od;
				
				if holes_filled then
					next;
				fi;
			else	# we have found something w\o the generating funciton
				x := res[1];
				y := res[2];
				checkholes := false;	# SkipY is good
				if extend_index <= 0 then	# update extend_index
					extend_index := nops(PPos1[i]) + 1;
				fi;
			fi;							

			# update the list
			PPos1[i] := [op(PPos1[i]), [x, y]];
						
			# update variables
			diffY := [op(diffY), y];
			diff_golden := [op(diff_golden), i + x - floor(y * golden)];
			X := [op(X), x, x + y];
#print(i, k1, [x, y], i + x - floor(y * golden));

			# Expand previous results if necessary
			res := TestExpansion(PPos1, i, x, y, X, SkipX, EXPANSION_SIZE);
			for l from 0 to i - 1 do
				PPos1[l] := res[1][l];
			od;
			if extend_index <= 0 then
				f1 := res[2];
				g := g - mtaylor(f1, [x1, x2], deg);
				f_extend := 0;
			else
				f_extend := f_extend + res[2];
			fi;
			X := res[3];
			SkipX := res[4];

			# test conditions for possible pattern
			
			# each y should be 1 plus the previous one
			if y = prev_y + 1 then
				if indx_y < 0 then
					indx_y := nops(PPos1[i]) - 1;
				fi;
			else
				indx_y := -10;
				prev_y := -10;
				checkholes := true;
			fi;
			prev_y := y;

			if checkholes or indx_y < 0 then
				next;
			fi;
			
			if x <= PPos1[i][indx_y][1] + PPos1[i][indx_y][2] then
				next;
			fi;
			
			# we want to test a number of consecutive pairs so that 
			# the second number of the first pair is greater than the first of the last pair
			for l from indx_y to nops(PPos1[i]) do
				if PPos1[i][l][1] + PPos1[i][l][2] > x then
					break;
				fi;
			od;
						
			if l <= nops(PPos1[i]) and l > MAX_EXT_SIZE and l >= pattern_start then
				res := TestPattern(diff_golden, l - MAX_EXT_SIZE, PPos1, i, 
					PatternIndx1, DiffY1, MIN_DIFFY_DIFF);
				
				if res < 0 then
					SkipX1[i] := SkipX;
					SkipY1[i] := SkipY;
					temp := diff_golden[l - MAX_EXT_SIZE..nops(PPos1[i])];
					DiffY1[i] := floor((max(op(temp)) + min(op(temp))) / 2);
					
					PatternIndx1[i] := [indx_y, -res, nops(PPos1[i])];
					
print(i, ` is good at `, PPos1[i], SkipX, SkipY, DiffY1[i], PatternIndx1[i]);
					break;
				else
					pattern_start := res;
if pattern_start > 0 then
print(`next pattern test start at`, pattern_start, PPos1[i][pattern_start]);
fi;
				fi;
			fi;
		od;
	od;
	
	PPos1, SkipX1, SkipY1, DiffY1, PatternIndx1;
end:

# recalculate generating function, due to new x, y, and deg
RecalculateGenFunc := proc(PPos, index, x, mex_Y, diffY, deg,
		MIN_CONVERT_EXT, checkholes, x1, x2)
	local time1, time2, f, f1, g1, g2, res, l;
	
	time1 := time();
	f := 1 / (1 - x1) / (1 - x2);
	f := f - ConvertAllPrevPPos(PPos, index, x, mex_Y, deg, x1, x2);
	res := ConvertSpecialPrevPPos(PPos, index, x1, x2);
	f := f - res[1];
	f1 := 0;
	
	for l from 1 to nops(PPos[index]) do
		if PPos[index][l][1] + PPos[index][l][2] < x
				and PPos[index][l][2] < mex_Y then	# what if the first < index????
			next;
		fi;

		f1 := f1 - ConvertOneCurrentPPos(PPos[index][l][1], PPos[index][l][2], x, mex_Y, deg, x1, x2);
	od;
					
	time2 := time();
	g2 := mtaylor(f, [x1, x2], deg);
	g1 := mtaylor(f1, [x1, x2], deg) + g2;
print(`Recalc generating function`, x, mex_Y, max(op(diffY)), `time = `, time2-time1, time()-time2, nops([op(f)]), deg, `time = `, time());

	g1, g2;
end:

# test if any of the previous results shall be extended
# basically x = mex(x, x+y), y = y + 1
TestExpansion := proc(PPos, index, x, y, X, SkipX, EXPANSION_SIZE)
	local i, l, tempPPos, res, f1, X1, SkipX1;
	
	i := index;
	
	for l from 0 to i - 1 do
		# if the second pile of the new set is getting close to 
		# that of the previous sets, expand the previous result sets.
		if x + i >= l + PPos[l][nops(PPos[l])][1] - 2 then
			tempPPos[l] := ExpandResult(PPos, l, x + EXPANSION_SIZE);
		else
			tempPPos[l] := [];
		fi;
	od;

	f1 := ConvertAllPrevPPos(tempPPos, i, -1, -1, x + y + EXPANSION_SIZE, x1, x2);
	res := ConvertSpecialPrevPPos(tempPPos, i, x1, x2);
	f1 := f1 + res[1];
	
	for l from 0 to i - 1 do
		tempPPos[l] := [op(PPos[l]), op(tempPPos[l])];
	od;
	X1 := [op(X), op(res[2])];
	SkipX1 := [op(SkipX), op(res[2])];

	tempPPos, f1, X1, SkipX1;
end:

# test if all the number not in diffY and less than y 
# are going to be permenantly skipped as the difference of X and Y
# return two items
# the first is true if all such number can be skipped, false otherwise
# the second is the list of number that will be skipped
TestSkippedNum := proc(index, x, y, diffY, f, deg, i_period, x1, x2)
	local i, j, g2, mex, tempdiffY, SkipY;
	
	SkipY := [];
	
	tempdiffY := diffY;
	mex := Mex(tempdiffY);
	
	while mex < y do
		
		g2 := coeff(f, x2, mex);
		# if there is no positive coeff starting form the current x 
		# to the specified deg, then the mex is skipped
		for i from min(x, deg - i_period - 2 - mex) to deg do
			if coeff(g2, x1, i) > 0 then
				RETURN(false, SkipY);
			fi;
		od;
		
		SkipY := [op(SkipY), mex];
		tempdiffY := [op(tempdiffY), mex];
		
		mex := Mex(tempdiffY);
	od;
	
	RETURN(true, SkipY);
end:

# Test if X list consist of numbers whose diffferences are at most 2
# and there are no consecutive numbers that are not the mean.
#
# return < 0 if it is good list, and the abs value is the start of the good pattern
# return 0 if the max difference > 2
# return the index for the end of the first set bad numbers if max difference <= 2
TestPattern := proc(difflist, start, PPos, index, PatternIndx1, DiffY1, MIN_DIFFY_DIFF)
	local i, j, l, maxd, mind, midd, count, newPPos, size, templist, tempdifflist,
		last_count, MAX_REPETENCE;

	MAX_REPETENCE := 3;
		
	size := nops(PPos[index]);

	# current P-position must have passed the start of the previous patterns
	for l from 0 to index - 1 do
#		if PPos[index][start][1] + index < PPos[l][PatternIndx1[l][2]][1] + l or
#				PPos[index][size][1] + index < PPos[l][PatternIndx1[l][3]][1] + l then
		if PPos[index][size][2] <= PPos[l][PatternIndx1[l][3]][2] then
			RETURN(0);
		fi;
	od;

	templist := difflist[start..size];
	maxd := max(op(templist));
	mind := min(op(templist));

	# the differences between the second piles and the golden section sequence
	# must be within three numbers
	if maxd - mind <= 2 then
		midd := floor((maxd + mind) / 2);
		count := 0;
		
		# the difference between this diffY and previous diffY 
		# must be at least MIN_DIFFY_DIFF
#		for l from 0 to index - 1 do
#			if DiffY1[l] - midd < MIN_DIFFY_DIFF then
			if index > 0 and DiffY1[index - 1] - midd < MIN_DIFFY_DIFF then
				RETURN(0);
			fi;
#		od;
		
		
		# no MAX_REPETENCE number of consecutive midd - 1 or midd + 1
		# if there are, we want to find the last occurrance
		last_count := nops(templist) + 1;
		l := 1;
		
		while l <= nops(templist) do
			for l from l to nops(templist) do
				if templist[l] <> midd then
					count := count + 1;
									
					if count = MAX_REPETENCE then
						break;
					fi;
				else
					count := 0;
				fi;
			od;
			
			count := 0;
			if l < nops(templist) then
				last_count := l;
			fi;
		od;
		l := last_count;
			
		if l > nops(templist) then
			# everything is good, trace back to the beginning of the pattern 
			# and stop here
			for l from start - 1 by -1 to 1 do
				if difflist[l] < midd - 1 or difflist[l] > midd + 1 then
					break;
				fi;
			od;
			l := l + 1;

			RETURN(-l);
		else
			for l from l to nops(templist) do
				if templist[l] = midd then
					break;
				fi;
			od;

			RETURN(start + l - 1);
		fi;
	fi;
	
	0;
end:

#############################################################################
# Find the next P-position without using the generating function
# the algorithm is proved in the associated paper
#############################################################################

FindNextPPos := proc(PPos, index, mex_X, diffY, mex_Y, DiffY1, PatternIndx1, MIN_DIFFY_DIFF)
	local i, j, golden, newdiff, size;
	
	golden := (1 + sqrt(5)) / 2:
	size := nops(PPos[index]);
	
	# the next x must fall into the patterns in previous results
	for i from 0 to index - 1 do
#		if mex_X < PPos[i][PatternIndx1[i][2]][1] then
		if mex_Y < PPos[i][PatternIndx1[i][2]][2] then
			RETURN(NULL);
		fi;
	od;
	
	# there can be no holes in diffY
	if  mex_Y <= max(op(diffY)) then
		RETURN(NULL);
	fi;

	# the difference of x and golden section pattern 
	# must be of some distance to the previous DiffY
	newdiff := index + mex_X - floor(mex_Y * golden):
	if index >= 1 then
		if DiffY1[index - 1] - newdiff < MIN_DIFFY_DIFF - 1 then
			RETURN(NULL);
		elif DiffY1[index - 1] - newdiff = MIN_DIFFY_DIFF - 1 then

			for i from 1 to nops(PPos[index-1]) do
				if PPos[index-1][i][2] = mex_Y - 1 then
					if index - 1 + PPos[index-1][i][1] <= index + mex_X then
						if index - 1 + PPos[index-1][i][1] - floor(golden * (mex_Y - 1)) 
							- DiffY1[index - 1] = -1 then
print([mex_X, mex_Y], `bad`, PPos[index-1][i][1] - floor(golden * (mex_Y - 1)) - DiffY1[index - 1],
	newdiff, DiffY1[index - 1]);
							RETURN(NULL);
						fi;
					fi;
					
					break;
				fi;
			od;			
#print([mex_X, mex_Y], `good`, PPos[index-1][i][1] - floor(golden * (mex_Y - 1)) - DiffY1[index - 1],
#	newdiff, DiffY1[index - 1]);
		fi;
	fi;

	RETURN([mex_X, mex_Y]);
end:

# expand previous results, assuming that the results are already "good",
# i.e. the second item of a pair in the list is one plus that of the previous pair,
# and it will keep it that way.
ExpandResult := proc(PPos, index, expandto)
	local i, j, X, diff, x, res, size;

	res := [];
	size := nops(PPos[index]);
	
	# make a list of all previously used numbers, inlcuding the skipped ones
	X := [seq(i, i = 0..PPos[index][size][1])];
	for i from 1 to size do
		X := [op(X), PPos[index][i][1] + PPos[index][i][2]];
	od;
	X := convert(convert(X, set), list);
	
	diff := PPos[index][size][2] + 1;
	x := Mex(X);
	
	# use the definition of X = mex(X, Y), and Y = X + n
	while x < expandto do
		res := [op(res), [x, diff]];
		
		X := [op(X), x, x + diff];
		x := Mex(X);
		diff := diff + 1;
	od;
	
	res;
end:

# find the nim sum (XOR) of two positive integers
NimSum := proc(x, y)
	local i, x1, y1, count, ret;

	x1 := x: 
	y1 := y: 
	count := 1: 
	ret := 0:
	
	while x1 * y1 > 0 do
		if (x1 mod 2) <> (y1 mod 2) then
			ret := ret + count;
		fi;

		count := count *2;
		x1 := floor(x1 / 2);
		y1 := floor(y1 / 2);
	od;

	ret + count * (x1 + y1);
end:

# find mex of a list of non-negtive numbers
Mex := proc(input)
	local i, inp, nstart, nend, nmid;
	
	inp := convert(convert(input, set), list);
	inp := sort(inp);
	
	nstart := 1;
	nend := nops(inp);
	
	if nend = 0 then
		RETURN(0);
	elif inp[nstart] > 0 then
		RETURN(0);
	elif inp[nend] = nend - 1 then
		RETURN(nend);
	fi;	

	while nend > nstart + 1 do 
		nmid := floor((nend + nstart) / 2);
		
		if nmid - 1 < inp[nmid] then
			nend := nmid;
		else 
			nstart := nmid;
		fi;		
	od;
	
	RETURN(nstart);
end: