# restart; read `\c:/Users/Asus/Google Drive/Aek/
# BreakGame/PartitionGame.txt`;

# First version: Jan 10, 2021

# This program accompanied my 
# own paper with Paul Ellis which 
# answers question 1 regarding to
# the nim sequence when C = {1,2c}, 
# c >=2 from the paper: 
# Partition Games. 


with(combinat,choose):
with(combinat,permute):


#################################
# Section 1,2: 
# Computational functions
#
# Mex(S), NimSum(S),
# Part(n,k), SubPart(n,k,L),
# NimVal(P,C), PGame(n,C),
# 
#################################

# Input: The list of non-negative
# integers S
# Output: Mex function of S, i.e.
# the first non-negative integer
# missing
# Try: Mex([0,1,3,6,2,4,12]);
Mex := proc(S) option remember;
local i,T;

T := sort([op(S)]);
i := 0;

while i+1 <= nops(T) and T[i+1] = i do
	i := i+1;
od:

i;	
end:	


# Input: The list of 
# non-negative integers S
# Output: NimSum of numbers in S
# Try: NimSum([7,5,3]);
NimSum := proc(S) option remember;
local i,s,T,t,m;

T := [seq(convert(s,base,2), s in S)];
m := max(seq( nops(t), t in T));

T := [seq( [op(t),0$(m-nops(t))] , t in T)];
T := [seq(add( t[i], t in T) mod 2,i=1..m)];

add( T[i]*2^(i-1),i=1..nops(T));
end:


# Input: non-negative integer n
# and positive integer p
# Output: Set of partitions of n
# with exactly p parts
# Try: Part(12,5);
Part := proc(n,p) option remember;
SubPart(n,p,n);
end:

# Input: non-negative integer n,
# positive integer p, 
# non-negative integer L
# Output: Set of partitions of n
# with exactly p parts,
# the largest part is no more than L
SubPart := proc(n,p,L) option remember;
local i,PP,T;

if p=0 then
	if n =0 then
		return({[]});
	else
		return({});
	fi:	
fi:

PP := {};
for i from 1 to L do
	PP := PP union 
	{seq([i,op(T)],T in SubPart(n-i,p-1,i))};
od:	

PP;	
end:


# Input: Partition P, a set 
# of positive integers C
# Output: Nim values of partition P
# according to C
# Try: NimVal([5,3,2],{1,4});
NimVal := proc(P,C) option remember;
local p;

NimSum([seq(PGame(p,C), p in P)]);
end:


# Input: non-negative integer n
# and a set of positive integers C
# Output: Grundy value of the heap 
# of size n with the cut C
# Try: [seq(PGame(n,{1,2,5}),n=1..36)];
# [seq(PGame(n,{1,6}),n=1..36)];
PGame := proc(n,C) option remember;
local c,P,p,G;

if member(0,C) then
	ERROR("Member of C must be positive");
fi:

if n <= min(op(C)) then 
	return(0);
fi:

G := {seq(seq(NimVal(P,C)
, P in Part(n,c+1)),c in C)};
Mex(G);
end:



#################################
# Section 3: 
# Some Simple Observations 
# and Nim-set
#
# CheckOB1(n,c), CheckOB2(c), 
# CheckOB3(n,a,c),
#
# NimSet(n,p,C), 
# Imp1(n,C), Imp2(n,C),
#
#################################

# Check my observation 1
# Input: positive integer n and c
# Output: Check observation 1
# Try:
# {seq(CheckOB1(n,3),n=1..50)}; 
# {seq(CheckOB1(n,5),n=1..150)};
CheckOB1 := proc(n,c) option remember;
local a,b;

if 0 = n mod 2*c or 4*c+1 = n mod 12*c then
	"Not Consider";
else
	a := CPGame(n+1,c);
	b := NimSum([CPGame(n,c),1]);
	evalb(a=b);
fi:	
end:


# Check my observation 2
# Input: positive integer n and c
# Output: Check observation 2
# Try:
# {seq(CheckOB2(c),c=2..40)}; 
CheckOB2 := proc(c) option remember;

[CNimVal([2*c,12*c],c),
CNimVal([4*c,10*c],c),
CNimVal([6*c,8*c],c),6];
	
end:


# Check my observation 3.2,3.3 and 3.4
# Input: positive integer n and c
# Output: Check observation 3.2,3.3 and 3.4
# Try:
# c:=2:{seq(seq(CheckOB3(n,a,c),a=1..40),n=1..20)}; 
CheckOB3 := proc(n,a,c) option remember;
local A,B;

A := [CNimVal([n,n,2*a*c],c),
CNimVal([n+1,n+1,2*a*c-2],c)];
A := evalb(A[1]=A[2]);

if 2 <> a mod 6 then
	B := [CNimVal([n+1,n+1,2*a*c+1],c),
	CNimVal([n,n,2*a*c+3],c)];
else
	B := [CNimVal([n+1,n+1,2*a*c+2],c),
	CNimVal([n,n,2*a*c+4],c)];	
fi:

B := evalb(B[1]=B[2]);
[A,B];
end:

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

# Input: non-negative integer n,
# positive integer p, 
# set of positive integers C
# Output: NimSet for breaking n 
# in to p piles and calculate
# nim value according to C 
# Try: 
# [seq(NimSet(n,2,{1,6}),n=1..30)];
# [seq(NimSet(n,7,{1,6}),n=1..30)];
# [seq(NimSet(n,3,{1,6}),n=1..26)];
# [seq(NimSet(n,4,{1,4,6}),n=1..30)];
# [seq(NimSet(n,9,{1,8}),n=9..40)];
NimSet := proc(n,p,C) option remember;
local P;
     
{seq(NimVal(P,C), P in Part(n,p))};
end:


# Important example 1:
# The value of the game is 
# the same as Nimset of 1 pile.
# Input: positive integer n
# and Cutset C
# Try: {seq(Imp1(n,{1,6}),n=1..30)};
Imp1 := proc(n,C) option remember;
local L,R;

L := {PGame(n,C)};
R := NimSet(n,1,C);
#print(L);
evalb(L=R);
end:


# Important example 2:
# The value of the game is Mex 
# of union of the nimsets.
# Input: positive integer n
# and Cutset C
# Try: {seq(Imp2(n,{1,2}),n=1..30)};
Imp2 := proc(n,C) option remember;
local c,L,R;

L := PGame(n,C);
R := Mex(`union`(seq(NimSet(n,c+1,C),c in C)));

#print(L,R);
evalb(L=R);
end:


#################################
# Section 4: 
# Development of Nimset
# 
# EntPar(n,p,c), 
# SubEntPart(n,p,L,c),
# ExitP(P,c), EntExit(N,p,c),
# GenNimSet(n,p,c),
#
#################################


# Entering partition! 
# Input: positive integers n,p,c
# Output: Entering partition 
# of n with exactly p parts
# consider the cut set {1,2c}
# Try: 
# EntPart(22,5,2);
# [seq(EntPart(n,2,3),n=1..30)];
EntPart := proc(n,p,c) option remember;
SubEntPart(n,p,n,c);
end:


# Input: positive integers n,p,L,c
# Output: Entering partition 
# of n with exactly p parts
# the largest part is no more than L
# consider the cut set {1,2c}
SubEntPart := proc(n,p,L,c) option remember;
local i,a,PP,T,S;

if p=0 then
	if n =0 then
		return({[]});
	else
		return({});
	fi:	
fi:

PP := {};
S := {seq(2*c*i+1,i=0..(L-1)/(2*c))}
union {seq(4*c+2+12*c*i,i=0..(L-4*c-2)/(12*c))};

for a in S do
	PP := PP union 
	{seq([a,op(T)],T in SubEntPart(n-a,p-1,a,c))};
od:	

PP;	
end:


# Input: parition P and 
# positive integer c
# Output: the exiting partition 
# of P consider the cut set {1,2c}
# Try: ExitP([1,1],3);
ExitP := proc(P,c) option remember;
local p,S;

S := [];
for p in P do
	if p mod 12*c = 4*c+1 then
		S := [op(S),p];
	elif p mod 12*c = 4*c+2 then
		S := [op(S),p+2*c-2];
	elif p mod 2*c = 1 then
		S := [op(S),p+2*c-1];
	else
		ERROR("p not supposed to be there");
	fi:
od:

S;	
end:


# Input: positive integers N, p and c
# Output: print the entering and
# the corresponding exiting partitions
# from 1 to N
# Try:
# EntExit(30,2,3);
EntExit := proc(N,p,c) option remember;
local n,P;

for n from 1 to N do
	for P in EntPart(n,p,c) do
		print(P,ExitP(P,c));
	od:
od:
return();
end:


# Generate the Nim set from the
# entering-exiting scheme.
# Implement the idea of Lemma 7,
# Corollary 9.
# Input: positive integers n, p and c
# Output: NimSet by generating from
# Entering partition
# Try: GenNimSet(30,2,3);
# p:=3: A:=[seq(NimSet(n,p,{1,8}),n=1..50)]:
# B:=GenNimSet(50,p,4); evalb(A=B);
GenNimSet := proc(n,p,c) option remember;
local i,t,P,T,Nim,A,S,K;

S := {};
Nim := [];

for i from 1 to n do
	for T in S do
		if T[3] < i then
			S := S minus {T};
		fi:
	od:	

	A := {seq([CNimVal(P,c),i,
		 add(t, t in ExitP(P,c))],
		 P in EntPart(i,p,c))};
	S := S union A; 
	K := {seq(NimSum([T[1],(i-T[2]) mod 2]), T in S)};
	Nim := [op(Nim),K];
od:	

return(Nim);
end:



#################################
# Section 4.1:
# Proof of Lemma 10: 
# Exiting partition
#
# ExitPart(n,p,c),
# SubExitPart(n,p,L,c),
# L10TestAll(c), 
#
# L10Proof(c), L10Case(P,c),
# L10Case1(P,c), L10Case2(P,c),
# L10Case3(P,c), L10Case4(P,c),
# Count1(P,h), FindMod(r,P,c);
#
#################################

# Exiting partition 
# Input: positive integers n, p and c
# Output: Exiting partition of n
# with exactly p parts consider
# cut-set {1,2c}  
# Try: ExitPart(142,5,3);
ExitPart := proc(n,p,c) option remember;
SubExitPart(n,p,n,c);
end:


# Input: positive integers n, p, L and c
# Output: Exit Partition of n
# with exactly p parts,
# the largest part is no more than L
# consider cut-set {1,2c}  
SubExitPart := proc(n,p,L,c) option remember;
local i,a,PP,T,S;

if p=0 then
	if n =0 then
		return({[]});
	else
		return({});
	fi:	
fi:

PP := {};
S := {seq(2*c*i,i=1..L/(2*c))}
union {seq(4*c+1+12*c*i,i=0..(L-4*c-1)/(12*c))};

for a in S do
	PP := PP union 
	{seq([a,op(T)],T in SubExitPart(n-a,p-1,a,c))};
od:	

PP;	
end:


# Check all the exiting cases 
# of lemma 10 (p=4). That is all 
# of the 7^4 cases of C={1,2c}
# Input: positive integer c >=2
# Output: P if there is no
# intermediate replacemnt for P
# and "Done" otherwise.
# Try:
# L10TestAll(2);
# L10TestAll(3);
L10TestAll := proc(c) option remember;
local i,S,h1,h2,h3,h4,T;

S := {seq(2*c*i,i=1..6),4*c+1};

for h1 in S do
for h2 in S do
for h3 in S do
for h4 in S do
	T := CMatchP([h1,h2,h3,h4],c);
	if T minus ExitPart(h1+h2+h3+h4,4,c) = {} then
		print([h1,h2,h3,h4]);
	fi:	
od: od: od: od:
return("Done");
end:


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


# Check the 4-case proof 
# of lemma 10 (p=4), C={1,2c}.
# Input: positive integer c >=2
# Output: P if there is no
# intermediate replacemnt for P
# and "Done" otherwise.
# Try:
# seq(L10Proof(c),c=2..7);
L10Proof := proc(c) option remember;
local i,S,h;

S := {seq(2*c*i,i=1..6),4*c+1};

for h[1] in S do
for h[2] in S do
for h[3] in S do
for h[4] in S do
	L10Case([h[1],h[2],h[3],h[4]],c);
od: od: od: od:	
return("Done");
end:


# Check for the intermediate 
# partition corresponding to P
# Input: partition P (4 parts)
# and positive integer c >=2
# Output: NULL if there is no 
# problem, otherwise print P.
L10Case := proc(P,c) option remember;
local n,p,q,case,Q;

n := add(p,p in P);
if Count1(P,4*c+1) =0 then
	Q := L10Case1(P,c);
	case := 1;
elif Count1(P,4*c+1) =1 then
	Q := L10Case2(P,c);
	case := 2;
elif Count1(P,4*c+1) =2 or Count1(P,4*c+1) =3 then
	Q := L10Case3(P,c);
	case := 3;
elif Count1(P,4*c+1) =4 then
	Q := L10Case4(P,c);
	case := 4;
else
	ERROR("NotSupposeLikeThis");
fi:
	
if add(q,q in Q) <> n or
CNimVal(Q,c) <> CNimVal(P,c) or
member(Q,ExitPart(n,nops(P),c))  then
	print("problem case",case, P,Q);
fi:

return();
end:


# Test case 1 of Lemma 10
# Input: partition P (4 parts)
# and positive integer c >= 2
# Output: non-exiting partition 
# Q (4 parts) with the same 
# nim-value as P
# Try: L10Case1([8,12,16,20],2);
L10Case1 := proc(P,c) option remember;
local T;

T := sort(P);

if nops({op(T)}) < 4 then
	if T[1]=T[2] then
		return([T[1]+1,T[2]+1,T[3]-2,T[4]]);
	elif T[2]=T[3] then
		return([T[1]-2,T[2]+1,T[3]+1,T[4]]);
	elif T[3]=T[4] then
		return([T[1]-2,T[2],T[3]+1,T[4]+1]);
	else
		ERROR("NotExpectThis");
	fi:	
elif  `subset`({2*c,12*c},{op(T)}) then
	if member(4*c,{op(T)}) or member(10*c,{op(T)}) then
		T := subs({2*c=4*c,12*c=10*c},T);
	else 	
		T := subs({2*c=6*c,12*c=8*c},T);
	fi:	
	L10Case1(T,c);
elif  `subset`({4*c,10*c},{op(T)}) then
	if member(2*c,{op(T)}) or member(12*c,{op(T)}) then
		T := subs({4*c=2*c,10*c=12*c},T);
	else 	
		T := subs({4*c=6*c,10*c=8*c},T);
	fi:	
	L10Case1(T,c);	
elif  `subset`({6*c,8*c},{op(T)}) then
	if member(2*c,{op(T)}) or member(12*c,{op(T)}) then
		T := subs({6*c=2*c,8*c=12*c},T);
	else 	
		T := subs({6*c=4*c,8*c=10*c},T);
	fi:	
	L10Case1(T,c);
else
	ERROR("NotExpectThis");
fi:	
end:


# Test case 2 of Lemma 10
# Input: partition P (4 parts)
# and positive integer c >= 2
# Output: non-exiting partition 
# Q (4 parts) with the same 
# nim-value as P
# Try: L10Case2([13,6,12,12],3);
L10Case2 := proc(P,c) option remember;
local T;

T := sort(P);

if nops({op(T)}) < 4 then
	if T[1]=T[2] and T[3] mod 2*c=0 then
		return([T[1]+1,T[2]+1,T[3]-2,T[4]]);
	elif T[1]=T[2] and T[4] mod 2*c=0 then
		return([T[1]+1,T[2]+1,T[3],T[4]-2]);	
	elif T[2]=T[3] and T[1] mod 2*c=0 then
		return([T[1]-2,T[2]+1,T[3]+1,T[4]]);
	elif T[2]=T[3] and T[4] mod 2*c=0 then
		return([T[1],T[2]+1,T[3]+1,T[4]-2]);	
	elif T[3]=T[4] and T[1] mod 2*c=0 then
		return([T[1]-2,T[2],T[3]+1,T[4]+1]);
	elif T[3]=T[4] and T[2] mod 2*c=0 then
		return([T[1],T[2]-2,T[3]+1,T[4]+1]);	
	else
		ERROR("NotExpectThis");
	fi:	
else
	T := sort([op({op(T)} minus {4*c+1})]);
	if T[1] = 4*c then
		return([6*c+1,2*c,T[2],T[3]]);
	elif T[1] = 6*c then
		return([8*c+1,2*c,T[2],T[3]]);
	elif T[1] = 12*c then
		return([10*c,6*c+1,T[2],T[3]]);	
	elif T[2] = 4*c then
		return([6*c+1,2*c,T[1],T[3]]);
	elif T[2] = 6*c then
		return([8*c+1,2*c,T[1],T[3]]);
	elif T[2] = 12*c then
		return([10*c,6*c+1,T[1],T[3]]);
	elif T[3] = 4*c then
		return([6*c+1,2*c,T[1],T[2]]);
	elif T[3] = 6*c then
		return([8*c+1,2*c,T[1],T[2]]);
	elif T[3] = 12*c then
		return([10*c,6*c+1,T[1],T[2]]);	
	else
		return([6*c+1,6*c,6*c,6*c]);
	fi:	
fi:

end:


# Test case 3 of Lemma 10
# Input: partition P (4 parts)
# and positive integer c >= 2
# Output: non-exiting partition 
# Q (4 parts) with the same 
# nim-value as P
# Try: L10Case3([6,13,13,18],3);
L10Case3 := proc(P,c) option remember;
local T,A,B;

A := FindMod(0,P,c);
B := FindMod(1,P,2*c);

if nops(A)=1 and nops(B)=3 then
	[B[1]+1,B[2]+1,B[3],A[1]-2];
elif nops(A)=2 and nops(B)=2 then
	[B[1]+1,B[2]+1,A[1]-2,A[2]];
else
	ERROR("NotSupposeLikeThis");
fi:	
end:


# Test case 4 of Lemma 10
# Input: partition P (4 parts)
# and positive integer c
# Output: intermediate partition 
# Q (4 parts) with the same 
# nim-value as P
# Try: L10Case4([13,13,13,13],3);
L10Case4 := proc(P,c) option remember;
[P[1]+1,P[2]+1,P[3]-1,P[4]-1];
end:


# Count number of element h in
# partition P
# Input: List P and number h
# Output: multiplicities of h in P
# Try: Count1([12,13,13,13],13);
Count1 := proc(P,h) option remember;
local c,p;

c := 0;
for p in P do
	if p=h then
		c := c+1;
	fi:
od:
c;
end:


# Input: integer r, list P and
# integer c >= 2
# Output: List of numbers in P 
# congruence to r mod 2*c 
# Try: FindMod(0,[13,13,18,24],3);
FindMod := proc(r,P,c)  option remember;
local p,T;

T := [];

for p in P do
	if p mod 2*c = r then
		T := [op(T),p];
	fi:
od:
T;	
end:



#################################
# Section 4.2: 
# Proof of Lemma 12
#
# L12TestAll(c),
# L12Proof(c), L12Case(P,c), 
# L12Case1(T), L12Case2(T,c),
# L12Case3(T,c),
#
#################################

# Check all the cases of lemma 12
# p=4, all the 7^4 cases of C={1,2c}
# Input: positive integer c >=2
# Output: P if there is no
# intermediate replacemnt for P
# and "Done" otherwise.
# Try:
# L12TestAll(2);
# L12TestAll(3);
L12TestAll := proc(c) option remember;
local i,S,h1,h2,h3,h4,T;

S := {seq(2*c*i+1,i=1..6),4*c+2};

for h1 in S do
for h2 in S do
for h3 in S do
for h4 in S do
	T := CMatchP([h1,h2,h3,h4],c);
	if T minus EntPart(h1+h2+h3+h4,4,c) = {} then
		print(sort([h1,h2,h3,h4]));
	fi:	
od: od: od: od:
return("Done");
end:



# Check the 3-case proof 
# of lemma 12 where p=4, C={1,2c}.
# Input: positive integer c >=2
# Output: P if there is no
# intermediate replacemnt for P
# and "Done" otherwise.
# Try:
# seq(L12Proof(c),c=2..7);
L12Proof := proc(c) option remember;
local a,S,h;

S := {seq(2*a*c+1,a=1..11),4*c+2};

for h[1] in S do
for h[2] in S do
for h[3] in S do
for h[4] in S do
	L12Case([h[1],h[2],h[3],h[4]],c);
od: od: od: od:	
return("Done");
end:


# Check for the non-entering
# partition corresponding to P
# Input: partition P (4 parts)
# and positive integer c >= 2
# Output: NULL if there is no 
# problem, otherwise print P.
L12Case := proc(P,c) option remember;
local i,n,q,loop,case,Q,T;

n := add(p,p in P);
loop := add(floor(P[i]/(12*c)),i=1..nops(P));
T := [seq(P[i] mod(12*c),i=1..nops(P))];
T := sort(T,`>`);

if T[1]=T[2] and T[3]=T[4] then
	Q := L12Case1(T,c);
	case := 1;
elif T[1]=T[2] or T[2]=T[3] or T[3]=T[4] then
	Q := L12Case2(T,c);
	case := 2;
elif T[1]<>T[2] or T[2]<>T[3] or T[3]<>T[4] then
	Q := L12Case3(T,c);
	case := 3;
else
	ERROR("L12Case");
fi:

Q := sort(Q,`>`);
i := floor(loop/4);
Q[1] := Q[1]+12*c*(i+BiAB(loop mod 4,3)); 
Q[2] := Q[2]+12*c*(i+BiAB(loop mod 4,2));
Q[3] := Q[3]+12*c*(i+BiAB(loop mod 4,1));
Q[4] := Q[4]+12*c*(i+BiAB(loop mod 4,0));
Q := sort(Q,`>`);

if add(q,q in Q) <> n or
CNimVal(Q,c) <> CNimVal(P,c) or 
member(0,{op(Q)}) or
member(Q,EntPart(n,nops(P),c))  then
	print("problem case",case, P,Q);	
fi:

return();
end:


# Test case 1 of Lemma 12
# Input: partition P (4 parts)
# and positive integer c >=2
# Output: non-entering partition 
# Q (4 parts) with the same 
# nim-value as P
# Try: L12Case1([31,31,11,11],5);
L12Case1 := proc(T,c) option remember;
[T[1]+1,T[2]+1,T[3]-1,T[4]-1];
end:


# Test case 2 of lemma 12
# Input: partition P (4 parts)
# and positive integer c >=2
# Output: non-entering partition 
# Q (4 parts) with the same 
# nim-value as P
# Try: L12Case2([7,7,13,14],3);
L12Case2 := proc(T,c) option remember;

if T[1]=T[2] and T[3] <> 4*c+1 then
	return([T[1]-1,T[2]-1,T[3]+2,T[4]]);
elif T[1]=T[2] and T[4] <> 4*c+1 then
	return([T[1]-1,T[2]-1,T[3],T[4]+2]);	
elif T[2]=T[3] and T[1] <> 4*c+1 then
	return([T[1]+2,T[2]-1,T[3]-1,T[4]]);
elif T[2]=T[3] and T[4] <> 4*c+1 then
	return([T[1],T[2]-1,T[3]-1,T[4]+2]);
elif T[3]=T[4] and T[1] <> 4*c+1 then
	return([T[1]+2,T[2],T[3]-1,T[4]-1]);
elif T[3]=T[4] and T[2] <> 4*c+1 then
	return([T[1],T[2]+2,T[3]-1,T[4]-1]);	
else
	ERROR("NotExpectThis");
fi:	

end:



# Test case 3 of Lemma 12
# Input: partition P (4 parts)
# and positive integer c >=2
# Output: non-entering partition 
# Q (4 parts) with the same 
# nim-value as P
# Try: L12Case3([25,14,7,1],3);
L12Case3 := proc(P,c) option remember;
local T,M;

T := P;
M := {op(P)};

# Case 3.1
if member(4*c+2,M) = false then
# Case 3.1a)
	if `subset`({8*c+1,10*c+1},M) then
		if member(1,M) then
			T := subs({1=2*c+1,10*c+1=8*c+1},T);
		elif member(4*c+1,M) then
			T := subs({4*c+1=6*c+1,10*c+1=8*c+1},T);
		elif member(6*c+1,M) then
			T := subs({6*c+1=4*c+1,8*c+1=10*c+1},T);
		else
			ERROR("NoCase 3.1 a)",P);
		fi:		
# Case 3.1b)	
	elif `subset`({1,6*c+1},M) then
		T := subs({1=2*c+1,6*c+1=4*c+1},T); 
# Case 3.1c)		
	elif T =[10*c+1,4*c+1,2*c+1,1] then
		T := [8*c+1,4*c+1,2*c+1,2*c+1];
	elif T =[8*c+1,6*c+1,4*c+1,2*c+1] then
		T := [10*c+1,4*c+1,4*c+1,2*c+1]; 
	elif T =[10*c+1,6*c+1,4*c+1,2*c+1] then
		T := [8*c+1,6*c+1,6*c+1,2*c+1]; 
	elif T =[8*c+1,4*c+1,2*c+1,1] then
		return([6*c,4*c,2*c+2,2*c+2]); 				
	else
		ERROR("NoCase 3.1 b)",P);
	fi:	
	
	T := sort(T,`>`);
	if T[1]=T[2] and T[3]=T[4] then
		return(L12Case1(T,c));
	elif T[1]=T[2] or T[2]=T[3] or T[3]=T[4] then
		return(L12Case2(T,c));
	else
		ERROR("NoCase 3.1, second",P);
	fi:	
# Case 3.2	4*c+2 is a member of M
else
	if member(4*c+1,M) then
		return(subs({4*c+2=8*c+2,4*c+1=1},T));
	elif member(6*c+1,M) then
		return(subs({4*c+2=8*c+2,6*c+1=2*c+1},T));	
	elif T = [10*c+1,8*c+1,4*c+2,2*c+1] then
		return([10*c,10*c,4*c+2,3]);
	elif T = [8*c+1,4*c+2,2*c+1,1] then
		return([8*c,3*c+1,3*c+1,3]);
	elif T = [10*c+1,4*c+2,2*c+1,1] then
		return([10*c,6*c,4,1]);		
	elif T = [10*c+1,8*c+1,4*c+2,1] then
		return([10*c+1,8*c,4*c,4]);	
	else
		ERROR("NoCase 3.2",P);
	fi:
fi:

end:


#################################
# Section 5: Lemma 16,17
#
# Lem16(N,p,c), MapPhi(n,p,c),
# PMod(n,c), 
# Lem17C(N,p,c),
# Lem16C2(n,c),
#
#################################

# Test the statement of Lemma 16
# Map from any c (>=3) to c=3.
# Input: positive integer N,p 
# and c (c>=3)
# Output: "Done" if 
# the nimset of C={1,2c},n=k+p-1 with p piles
# equals the nimset of C={1,6},n=phi_p(k+p-1) 
# with p piles from n=p..N.
# Otherwise print a problem case.
# Try: 
# Lem16(6,8,4);
# seq(Lem16(170,p,5),p=2..6);
Lem16 := proc(N,p,c) option remember;
local m,n,A,B;

if c < 3 then
	ERROR("c is too small");
fi:	

for n from p to N do
	A := GenNimSet(n,p,c)[n];
	m := MapPhi(n,p,c);
	B := GenNimSet(m,p,3)[m];

	if evalb(A=B) = false then
		print([n,A,m,B]);
	fi:	
od:	

return("done");
end:


# Mapping of phi from 
# definition 14
# Input: positive integer n,p 
# and c (c>=3)
# Output: phi_p(n+p-1) according
# to definition 14
# Try: 
# [seq(MapPhi(n,3,4),n=3..10)];
# [seq(MapPhi(n,1,4),n=1..80)];
MapPhi := proc(n,p,c) option remember;
local r,q,rp;

if c < 3 then
	ERROR("c is too small");
elif n = 1 then
	return(1);
elif n < p then
	ERROR("n is less than p");
fi:	

r := PMod(n-p+1,2*c);
q := (n-(r+p-1))/(2*c);

if r >= 2*c-1 then 
	rp := 6-(r mod 2);
elif r = 1 or r = 2 then 
	rp := r;	
else
	rp := 4-(r mod 2):
fi:	

6*q+rp+p-1;
end:

# Modulo 
# Input: integers n and c
# Output: n mod c with 
# k*c mod c = c (instead of 0).
# Try: [seq(PMod(n,8),n=-20..20)];
PMod := proc(n,c) option remember;
local r;

r := n mod c;
if r = 0 then
	c;
else
	r;
fi:
end:


# Lemma 17c):
# Test the linear mapping for 
# specific p=2,3,4 and c
# (with one exception of case p=4).
# Input: positive integers N,p 
# and c (c>=3)
# Output: "Done" if 
# all the linear map of 
# entering/exiting partitions 
# for fix p and c works, n=1..N.
# Otherwise print a problem case.
# Try: 
# Lem17C(500,2,5);
# Lem17C(500,3,5);
# Lem17C(500,4,5);
Lem17C := proc(N,p,c) option remember;
local n,t,P;

if c < 3 then
	ERROR("c is too small");
elif p<2 or p>4 then
	ERROR("p must be 2 or 3 or 4");
fi:

for n from 1 to N do
for P in EntPart(n,p,c) do
	if add(MapPhi(t,1,c),t in P) 
	<> MapPhi(n,p,c) then
		print(n,P);
	fi:
od:	

for P in ExitPart(n,p,c) do
	if add(MapPhi(t,1,c),t in P) 
	<> MapPhi(n,p,c) then
		print(n,P);
	fi:
od:	
od:

print("Done");
end:


# Check case p=4 of Lemma 16
# Input: integer n and c
# Output: "Done" if the claim
# for n and c in this case 
# is good. Otherwise print
# the problematic case.
# Try: 
# seq(Lem16C2(n,5),n=1..64);
# seq(Lem16C2(n,4),n=1..64);
Lem16C2 := proc(n,c) option remember;
local s,t,S,T;

for S in EntPart(n,4,c) do
if member(1,{op(S)}) then
	T := [seq(MapPhi(s,1,c), s in S)];
	if NimVal(S,{1,2*c}) <> NimVal(T,{1,6})
	or MapPhi(n,4,c) <> add(t, t in T) then
		return(S,T,NimVal(S,{1,2*c})
		,NimVal(S,{1,2*c}));
	fi:	
fi:	
od:
return("Done");
end:


################################
# Section 5.1: Theorem 15 
# 
# Thm15Ca(N,c), Thm15Cb(N,c),
#
#################################

# Check claim (a) of Theorem 15
# Input: integer N and c
# Output: "Done" if the claims
# for n=2..N and c are good. 
# Otherwise print the problematic case.
# Try: 
# Thm15Ca(230,5);
Thm15Ca := proc(N,c) option remember;
local m,n,A,B;

A := GenNimSet(N,2*c+1,c);
m := MapPhi(N,1,c);
B := GenNimSet(m,7,3);

for n from 2 to N do
	m := MapPhi(n,1,c);
	if evalb(A[n]=B[m]) = false then
		print(n,A[n],m,B[m]);
	fi:	
od:	

return("Done");
end:


# Check claim (b) of Theorem 15
# Input: integer N and c
# Output: "Done" if the claims
# for n=2..N and c are good. 
# Otherwise print the problematic case.
# Try: 
# Thm15Cb(330,5);
Thm15Cb := proc(N,c) option remember;
local m,n,A,B;

for n from 2 to N do
	A := GenNimSet(n,2,c)[n];
	m := MapPhi(n,1,c);
	B := GenNimSet(m,2,3)[m];

	if evalb(A=B) = false then
		print(n,A,m,B);
	fi:	
od:	

return("Done");
end:


################################
# Section 6: An extension
# 
# Lem18(n,p,c), Thm19(n,C),
# 
#################################

# Test Lemma 18
# Input: integers n,p and c
# Output: The set of nim-values
# of partition of n with p parts
# under the cut {1,2c}.
# Try:
# Lem18(10,3,3);
# seq(print(seq(Lem18(51,p,c),c=2..5)),p=4..6);
Lem18 := proc(n,p,c) option remember;
local P,L,R;

R := {seq(CNimVal(P,c),P in Part(n,p))};
L := {seq(CNimVal(P,c),P in Part(n,p+2))};

[L,L minus R];
end:


# Input: integer n and set C
# Output: The nim-values of n under 
# the cut {1,2*C} and {1,2*minC}.
# Try:
# Thm19(30,{3,5,7});
Thm19 := proc(N,C) option remember;
local a,c,n,L,R;

c := min(op(C));
if c < 2 then
	ERROR("Small c");
fi:	

for n from 1 to N do
	L := PGame(n,{1,seq(2*a,a in C)});
	R := PGame(n,{1,2*c});

	if L <> R then
		print("Wrong",L,R);
	fi:
od:

return("Done");
end:


################################
# Section 7: What Else?
# 
# Conj1(N,S), Conj2(n,S),
# 
#################################

# Input: integer N and set S
# Output: the list of values 
# of the game from 1 to N 
# according to set S restricted
# by condition of conjecture 1
# Try: Conj1(20,{4,5,6});
Conj1 := proc(N,S) option remember;
local n,s,X,Y;

X :={};
Y := {};

if min(op(S)) <=3 then
	return("Member must be at least 4");
fi:	
	
for s in S do
	if s mod 2 = 0 then
		X := X union {s};
	else
		Y := Y union {s};
	fi:
od:	

[seq(PGame(n,{1,3} union X union Y),n=1..N)];
end:


# Input: integer N and set S
# Output: the list of values 
# of the game from 1 to N 
# according to set S restricted
# by condition of conjecture 2
# Try: 
# Conj2(20,{4,8,13});
# Conj2(20,{6,19});
Conj2 := proc(N,S) option remember;
local n,s,X,Y,x,y,A,B;

X :={};
Y := {};

if min(op(S)) <=3 then
	return("Member must be at least 4");
fi:	
	
for s in S do
	if s mod 2 = 0 then
		X := X union {s};
	else
		Y := Y union {s};
	fi:
od:	

x := min(op(X));
y := min(op(Y));

if 3*x >= y then 
	return("NotTheCase"):
fi:

A:=[seq(PGame(n,{1} union X union Y),n=1..N)];
B:=[seq(PGame(n,{1,x}),n=1..N)];

return(A,evalb(A=B));
end:


#################################
# Section 8: Helper functions 
# 
# CPGame(n,c), CNimVal(P,c),
# CMatchP(P,c),
# BiAB(a,b),
#
#################################

# Fast computation
# Input: non-negative integer n,  
# a positive integer c >= 2
# Output: Value of PGame(n) 
# of C={1,2c} from the conjectures
# Try: [seq(CPGame(n,3),n=1..48)];
# A:=[seq(CPGame(n,3),n=1..48)];
# B:=[seq(PGame(n,{1,6}),n=1..48)]:
# evalb(A=B);
CPGame := proc(n,c) option remember;
local q,r,t,Ans;

if c =1 or c=0 or n < 1 then
	ERROR("NoCase",n,c);
fi:	

r := PMod(n,2*c);
q := (n-r)/(2*c);
t := q mod 6;
Ans := [[0,1],[2,3],[5,4],[3,2],[4,5],[6,7]];

if t=2 and r=1 then
	(q-t)/6*8+1;
else
	(q-t)/6*8+Ans[t+1,PMod(r,2)];
fi:
end:


# Fast computation
# Input: Partition P,  
# a positive integer c
# Output: Nim values of the 
# partition P with C={1,2c}
# Try:
# CNimVal([5,3,2],2);
# CNimVal([5,13,2],2);
# NimVal([5,13,2],{1,4});
CNimVal := proc(P,c) option remember;
local p;

NimSum([seq(CPGame(p,c), p in P)]);
end:


# Fast computation
# Input: the list P i.e. partition
# of n with p parts
# Output: the set of partitions of 
# n with p part with the same 
# nim-values as P.
# Try: CMatchP([7,7,7,7,7],6);
CMatchP := proc(P,c) option remember;
local p,T,n,A;

n := add(p, p in P);
A := {};

for T in Part(n,nops(P)) do
	if CNimVal(T,c) = CNimVal(P,c) then
		A := A union {T};
	fi:
od:
	
return(A);
end:


# Input: integers a and b
# Output: 1 if a>b and 0 otherwise
# Try: BiAB(5,2);
BiAB := proc(a,b) option remember;
if a > b then 
	return(1);
else 
	return(0);
fi:

end: