######################################################################
#JGgraph.txt Save this file as JGgraph.txt                           #
#To use it, stay in the                                              #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read JGgraph.txt : <Enter>                         #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger & Blair Seidler, Rutgers University ,  #
#doronzeil at gmail dot com                                          #
#blair at math dot rutgers dot edu                                   #
######################################################################
 
#Created: Nov. 28, 2022
 
print(`Created: Nov. 28, 2022`):
print(`Current version: Jun 23, 2023`):
print(``):
print(`This is JuniperGreen, a Maple package designed to analyze`):
print(`the Juniper Green game.`):
print(``):
print(`Procedures with "h" in their names refer to the "hard" version`):
print(`of the game in which the first player must choose an even number.`):
print(``):
print(`Procedures with "g" in their names refer to a variant with two`):
print(`additional parameters A and B. The legal moves are divisors of`):
print(`the current number n which are at least n/A and multiples of the`):
print(`current number which are at most n*B`):
print(``):
print(`Procedures with "a" in their names refer to a variant with two`):
print(`additional parameters A and B. The legal moves are found by`):
print(`subtracting an element of set A or adding an element of set B.`):
print(``):
print(`written by Doron Zeilberger and Blair Seidler.`):
print(``):
print(`Please report bugs to blair at math dot rutgers dot edu`):
print(``):
print(`For a list of the procedures type Help();, for help with`):
print(`a specific procedure, type Help(procedure_name);   .`):
print(``):

with(combinat):
with(numtheory):
with(GraphTheory):

Help:=proc():

if args=NULL then
print(`The main procedures are: `):
print(`M,WM,G,Mh,WMh,Gh,Mg,WMg,Gg,Mhg,WMhg,Ghg,Ma,WMa,Ga`):
print(``):

elif nops([args])=1 and op(1,[args])=M then
print(`M(n,P): All the legal moves in Juniper Green with starting integer n,`):
print(` P=[m,S], where m is current integer m and "already used set" S.`):
print(``):
print(` Try: M(9,[9,{}]);`):

elif nops([args])=1 and op(1,[args])=WM then
print(`WM(n,P): The set of winning moves in Juniper-Green starting with n,`):
print(` P=[m,S], where m is current integer m and "already used set" S.`):
print(``):
print(` Try: WM(9,[9,{}]);`):

elif nops([args])=1 and op(1,[args])=G then
print(`G(n,P): The Sprague-Grundy value of position P in EASY`):
print(`Juniper Green with starting integer n.`):
print(``):
print(` Try: G(9,[9,{}]);`):

elif nops([args])=1 and op(1,[args])=Mh then
print(`Mh(n,P): All the legal moves in Juniper Green with starting integer n,`):
print(` P=[m,S], where m is current integer m and "already used set" S.`):
print(``):
print(` Try: Mh(9,[9,{}]);`):

elif nops([args])=1 and op(1,[args])=WMh then
print(`WMh(n,P): The set of winning moves in Juniper-Green starting with n,`):
print(` P=[m,S], where m is current integer m and "already used set" S.`):
print(``):
print(` Try: WMh(9,[9,{}]);`):

elif nops([args])=1 and op(1,[args])=Gh then
print(`Gh(n,P): The Sprague-Grundy value of position P in HARD`):
print(` Juniper Green with starting integer n.`):
print(``):
print(` Try: Gh(9,[9,{}]);`):

elif nops([args])=1 and op(1,[args])=mex then
print(`mex(A): Returns least nonnegative integer not in set A.`):
print(``):
print(` Try: mex({0,1,3,4,6});`):

elif nops([args])=1 and op(1,[args])=Mg then
print(`Mg(n,P,A,B): All the legal moves in Juniper Green with starting integer n,`):
print(` P=[m,S], where m is current integer m and "already used set" S. In this`):
print(` variation, all legal moves must be between m/A and m*B`):
print(``):
print(` Try: Mg(9,[9,{}],2,2);`):

elif nops([args])=1 and op(1,[args])=WMg then
print(`WMg(n,P,A,B): The set of winning moves in Juniper-Green starting with n,`):
print(` P=[m,S], where m is current integer m and "already used set" S. In this`):
print(` variation, all legal moves must be between m/A and m*B`):
print(``):
print(` Try: WMg(9,[9,{}],2,2);`):

elif nops([args])=1 and op(1,[args])=Gg then
print(`Gg(n,P): The Sprague-Grundy value of position P in EASY`):
print(` Juniper Green with starting integer n. In this`):
print(` variation, all legal moves must be between m/A and m*B`):
print(``):
print(` Try: Gg(9,[9,{}],2,2);`):

elif nops([args])=1 and op(1,[args])=Mhg then
print(`Mhg(n,P,A,B): All the legal moves in Juniper Green with starting integer n,`):
print(` P=[m,S], where m is current integer m and "already used set" S. In this`):
print(` variation, all legal moves must be between m/A and m*B`):
print(``):
print(` Try: Mhg(9,[9,{}],2,2);`):

elif nops([args])=1 and op(1,[args])=WMhg then
print(`WMhg(n,P,A,B): The set of winning moves in Juniper-Green starting with n,`):
print(` P=[m,S], where m is current integer m and "already used set" S. In this`):
print(` variation, all legal moves must be between m/A and m*B`):
print(``):
print(` Try: WMhg(9,[9,{}],2,2);`):

elif nops([args])=1 and op(1,[args])=Ghg then
print(`Ghg(n,P): The Sprague-Grundy value of position P in HARD`):
print(` Juniper Green with starting integer n. In this`):
print(` variation, all legal moves must be between m/A and m*B`):
print(``):
print(` Try: Ghg(9,[9,{}],2,2);`):

elif nops([args])=1 and op(1,[args])=Ma then
print(`Ma(n,P,A,B): All the legal moves in additive Juniper Green with starting integer n,`):
print(` P=[m,S], where m is current integer m and "already used set" S. In this`):
print(` variation, all legal moves must be of the form m-a for some a in set A`):
print(` or form m+b for some b in set B.`):
print(``):
print(` Try: Ma(9,[4,{4}],{2},{2});`):

elif nops([args])=1 and op(1,[args])=WMa then
print(`WMa(n,P,A,B): The set of winning moves in additive Juniper Green starting with n,`):
print(` P=[m,S], where m is current integer m and "already used set" S. In this`):
print(` variation, all legal moves must be of the form m-a for some a in set A`):
print(` or form m+b for some b in set B.`):
print(``):
print(` Try: WMa(9,[9,{}],{2},{2});`):

elif nops([args])=1 and op(1,[args])=Ga then
print(`Ga(n,P): The Sprague-Grundy value of position P in additive`):
print(` Juniper Green with starting integer n. In this`):
print(` variation, all legal moves must be of the form m-a for some a in set A`):
print(` or form m+b for some b in set B.`):
print(``):
print(` Try: Ga(9,[9,{}],{2},{2});`):

else print(`No help for `,op(1,[args])):
fi:

end:



mex:=proc(A) local i: 
if A={} then
 RETURN(0):
fi:
for i from 0 while member(i,A) do od:i:end:


#EASY VERSION
#M(n,P): All the legal moves in Juniper Green with starting integer n, P=[m,S], where m is current integer m and "already used set" S. Try:
#M(10,[10,{}]);
M:=proc(n,P) local m,S,i,S1,m1:
option remember:
m:=P[1]:
S:=P[2]:

if P[1]=n and S={} then
S1:={seq(i,i=1..n)}:
else
S1:=divisors(m) union {seq(m*i,i=1..trunc(n/m))}:
fi:

S1:=S1 minus S:
{seq([m1, S union {m1}],m1 in S1)}:
end:

#WM(n,P): The set of winning moves in Juniper-Green starting with n and current state P. Try:
#WM(10,[10,{}]);
WM:=proc(n,P)  local gu,gu1,mu:
option remember:
gu:=M(n,P):
if gu={} then
 RETURN({}):
fi:

mu:={}:

for gu1 in gu do
 if WM(n,gu1)={} then
  mu:=mu union {gu1}:
fi:
od:
mu:


end:

#G(n,P): The Sprague-Grundy value of position P in EASY Juniper Green with starting integer n
G:=proc(n,P) local S,s:
option remember:

S:=M(n,P):

if S={} then
 0:
else
mex({seq(G(n,s),s in S)}):
fi:


end:

#HARD VERSION

#Mh(n,P): All the legal moves in Juniper Green Hard with starting integer n, P=[m,S], where m is current integer m and "already used set" S. Try:
#Mh(10,[10,{}]);
Mh:=proc(n,P) local m,S,i,S1,m1:
option remember:
m:=P[1]:
S:=P[2]:

if P[1]=n and S={} then
S1:={seq(2*i,i=1..trunc(n/2))}:
else
S1:=divisors(m) union {seq(m*i,i=1..trunc(n/m))}:
fi:

S1:=S1 minus S:
{seq([m1, S union {m1}],m1 in S1)}:
end:

(* Original version
Mhgraph:=proc(n) local G,i,j,v:
option remember:
G:=Digraph(n):
for v from 1 to n do
#print(G,{seq([v,i],i in divisors(v) union {seq(v*i,i=1..trunc(n/v))})}):
  AddArc(G,{seq([v,i],i in (divisors(v) union {seq(v*i,i=1..trunc(n/v))}) minus {v})}):
od:
G:
end:
*)

Mhgraph:=proc(n) local G,i,j,v:
option remember:
G:=table():
for v from 1 to n do
#print(G,{seq([v,i],i in divisors(v) union {seq(v*i,i=1..trunc(n/v))})}):
  G[v]:={seq(i,i in (divisors(v) union {seq(v*i,i=1..trunc(n/v))}) minus {v})}:
od:
op(G):
end:

Ghgraph1:=proc(n,c,S) local N,i:
option remember:
N:=Mhgraph(n)[c] minus S:
if nops(N)=0 then
 0:
else
 mex({seq(Ghgraph1(n,i,S union {i}),i in N)}):
fi:
end:

Ghgraph:=proc(n) local G,v,sv,i,j:
option remember:
sv:={seq(2*j,j in 1..trunc(n/2))}:
mex({seq(Ghgraph1(n,i,{i}),i in sv)}):
end:

#WMh(n,P): The set of winning moves in Juniper-Green starting with n and current state P. Try:
#WMh(10,[10,{}]);
WMh:=proc(n,P)  local gu,gu1,mu:
option remember:
gu:=Mh(n,P):
if gu={} then
 RETURN({}):
fi:

mu:={}:

for gu1 in gu do
 if WMh(n,gu1)={} then
  mu:=mu union {gu1}:
fi:
od:
mu:


end:

#Gh(n,P): The Sprague-Grundy value of position P in Hard Juniper Green with starting integer n
Gh:=proc(n,P) local S,s:
option remember:

S:=Mh(n,P):

if S={} then
 0:
else
mex({seq(Gh(n,s),s in S)}):
fi:


end:

#Variations with A, B
#This is a version where a legal move is
#m->m/a where a<=A and a is divisible by m
#or
#m->m*a where a<=B
#of  course not using any integer used before.


#EASY VERSION
#Mg(n,P,A,B): All the legal moves in Juniper Green with starting integer n, P=[m,S], where m is current integer m and "already used set" S. Try:
#Mg(10,[10,{}],2,2);
Mg:=proc(n,P,A,B) local m,S,i,S1,m1:
option remember:
m:=P[1]:
S:=P[2]:

if P[1]=n and S={} then
S1:={seq(i,i=1..n)}:
else
S1:=(divisors(m) intersect {seq(m/i,i=2..A)}) union {seq(m*i,i=1..min(B,trunc(n/m)))}:
fi:

S1:=S1 minus S:
{seq([m1, S union {m1}],m1 in S1)}:
end:

#WMg(n,P,A,B): The set of winning moves in Juniper-Green starting with n and current state P. Try:
#WMg(10,[10,{}],2,2);
WMg:=proc(n,P,A,B)  local gu,gu1,mu:
option remember:
gu:=Mg(n,P,A,B):
if gu={} then
 RETURN({}):
fi:

mu:={}:

for gu1 in gu do
 if WMg(n,gu1,A,B)={} then
  mu:=mu union {gu1}:
fi:
od:
mu:


end:

#Gg(n,P,A,B): The Sprague-Grundy value of position P in EASY Juniper Green with starting integer n
Gg:=proc(n,P,A,B) local S,s:
option remember:

S:=Mg(n,P,A,B):

if S={} then
 0:
else
mex({seq(Gg(n,s,A,B),s in S)}):
fi:


end:

#HARD VERSION
#Mhg(n,P,A,B): All the legal moves in Juniper Green with starting integer n, P=[m,S], where m is current integer m and "already used set" S. Try:
#Mhg(10,[10,{}],2,2);
Mhg:=proc(n,P,A,B) local m,S,i,S1,m1:
option remember:
m:=P[1]:
S:=P[2]:

if P[1]=n and S={} then
S1:={seq(2*i,i=1..trunc(n/2))}:
else
S1:=(divisors(m) intersect {seq(m/i,i=2..A)}) union {seq(m*i,i=1..min(B,trunc(n/m)))}:
fi:

S1:=S1 minus S:
{seq([m1, S union {m1}],m1 in S1)}:
end:

#WMhg(n,P,A,B): The set of winning moves in Juniper-Green starting with n and current state P. Try:
#WMhg(10,[10,{}],2,2);
WMhg:=proc(n,P,A,B)  local gu,gu1,mu:
option remember:
gu:=Mhg(n,P,A,B):
if gu={} then
 RETURN({}):
fi:

mu:={}:

for gu1 in gu do
 if WMhg(n,gu1,A,B)={} then
  mu:=mu union {gu1}:
fi:
od:
mu:


end:

#Ghg(n,P,A,B): The Sprague-Grundy value of position P in EASY Juniper Green with starting integer n
Ghg:=proc(n,P,A,B) local S,s:
option remember:

S:=Mhg(n,P,A,B):

if S={} then
 0:
else
mex({seq(Ghg(n,s,A,B),s in S)}):
fi:


end:


#Additive Variation with sets A, B
#This is a version where a legal move is
#m->m-a where a in A or
#m->m+b where b in B
#of  course not using any integer used before.

#Ma(n,P,A,B): All the legal moves in additive Juniper Green with starting integer n, P=[m,S], where m is current integer m and "already used set" S. Try:
#Ma(10,[10,{}],{1,2},{1,2});
Ma:=proc(n,P,A,B) local m,S,i,S1,m1:
option remember:
m:=P[1]:
S:=P[2]:

if P[1]=n and S={} then
S1:={seq(i,i=1..n)}:
else
S1:=({seq(m-m1,m1 in A)} union {seq(m+m1,m1 in B)}) intersect {seq(i,i=1..n)}:
fi:

S1:=S1 minus S:
{seq([m1, S union {m1}],m1 in S1)}:
end:

Magraph:=proc(n,A,B) local G,i,j,v,U:
option remember:
G:=table():
U:={seq(i,i=1..n)}:
for v from 1 to n do
  G[v]:={seq(i,i in ({seq(v-j,j in A)} union {seq(v+j,j in B)}) intersect (U minus {v}))}:
od:
op(G):
end:

Gagraph1:=proc(n,A,B,c,S) local N,i:
option remember:
N:=Magraph(n,A,B)[c] minus S:
if nops(N)=0 then
 0:
else
 mex({seq(Gagraph1(n,A,B,i,S union {i}),i in N)}):
fi:
end:

Gagraph:=proc(n,A,B) local G,v,sv,i,j:
option remember:
mex({seq(Gagraph1(n,A,B,i,{i}),i=1..n)}):
end:


#WMa(n,P,A,B): The set of winning moves in additive Juniper-Green starting with n and current state P. Try:
#WMa(10,[10,{}],2,2);
WMa:=proc(n,P,A,B)  local gu,gu1,mu:
option remember:
gu:=Ma(n,P,A,B):
if gu={} then
 RETURN({}):
fi:

mu:={}:

for gu1 in gu do
 if WMa(n,gu1,A,B)={} then
  mu:=mu union {gu1}:
fi:
od:
mu:


end:

#Ga(n,P,A,B): The Sprague-Grundy value of position P in EASY additive Juniper Green with starting integer n
Ga:=proc(n,P,A,B) local S,s:
option remember:

S:=Ma(n,P,A,B):

if S={} then
 0:
else
mex({seq(Ga(n,s,A,B),s in S)}):
fi:


end:


#4. Write a procedure FindUP(L) that inputs a list L and outputs a pair of lists L1,L2 such that
# (conjecturally) L=[op(L1),op(L2)^infinity]
# i.e., after an initial segment L1, it starts being periodic of period nops(L2)
FindUP:=proc(L) local i,n,T,P:
n:=nops(L):
for i from 1 to trunc(n/2) do
  T:=[op(i..n,L)]:
  P:=FindPer(T):
  if type(P,list) then RETURN([[op(1..i-1,L)],P]): fi:
od:
FAIL:
end:

#5. By using FindUP(L) and Ga(n,[n,{}],A,B), find(if possible) conjectured ultimate periodic descriptions
# of [seq(Ga(n,s,A,B),b=1..infinity)]
FindConjPerGa:=proc(S,n,step) local i,b,P,s,s1:
s:=[]:
for i from step to n by step do
  s1:=[seq(Ga(b,[b,{}],S,S), b = i-step+1..i)]:
  s:=[op(s),op(s1)]:
  P:=FindUP(s):
  if P<>FAIL then
    print(P):
    printf("set=%s n=%d, Initial segment: %d, Period: %d\n",convert(S,string),i,nops(P[1]),nops(P[2])):
    break:
  fi:
od:
end:

#IsPer1(L,p): Is the list L periodic of period p
IsPer1:=proc(L,p) local i:
for i from 1 to nops(L)-p do
 if L[i]<>L[i+p] then
  RETURN(false):
fi:
od:
true:
end:

#FindPer(L): finds the smallest period or return FAIL
FindPer:=proc(L) local p:
for p from 1 to trunc(nops(L)/2) do
if IsPer1(L,p) then
   RETURN([op(1..p,L)]):
fi:
od:
FAIL:
end: