#returns prob of winning.  if doesAhave9h = 0, then returns prob of winning 
#given A does not have 9h.  Otherwise, assume A has 9h. k must be <= 4.
probofwinning:= proc(k, doesAhave9h)

local i, S, numoccurs, j, winlose, prob, numwins, numdeals, index;

S:=startpos(k);

numoccurs:= [seq(binomial(3,S[i][1][1])*product(binomial(4,S[i][1][j]),j=2..k),i=1..nops(S))];
#numoccurs is the number of ways in which the given starting hand can be dealt

if k = 1 then
 winlose:= winorlose(1,4);
elif k = 2 then
 winlose:= winorlose(2,16);
elif k = 3 then
 winlose := winorlose(3,39);
elif k = 4 then
 winlose := winorlose(4,85);
else return("ERROR");
fi;

numwins:= 0;
numdeals:= 0;
for i from 1 to nops(S) do
 if (doesAhave9h = 0 and winlose[S[i]] = 1) or (doesAhave9h <> 0 and winlose[S[i]]=-1) then
  numwins:= numwins + numoccurs[i];
 fi;
 numdeals:=numdeals + numoccurs[i];
od;


return((numwins/numdeals));


end proc;



generatemoves:= proc(k)

local i, S, allchildren, pos;
S:= {op(startpos(k))};

allchildren:={};

do 
 for pos in S do
  allchildren:=allchildren union {op(legalmoves(op(pos)))};
 od;
 if allchildren subset S then 
  return(S);
 fi;
 S:= S union allchildren;
od;

end proc;


#Find a sequence from a starting position to (A,B) based on the claim and corollary
printpath:= proc(A,B) 
 local k, i,j, currpos, lowest, highest, temppile, l, t, q, r, isAsturn, numberofturns, n, a,b;
 k := nops(A);
 a := sum(A[i],i=1..k);
 b := sum(B[j],j = 1..k);
 if a < 1 or b < 3 then return("ERROR, player to move must have at least 1 card and other player must have at least 3 cards.") fi;
 n := 4*k - 1 - a - b;


if n = 0 then
 numberofturns :=0;
 isAsturn := 0;
 lowest:= 0;
 currpos := [[seq(2, i=1..k)], [1, seq(2, j =1..k-1)]];
 print(currpos);

 i := 1;
 while i <= k do   #keep going through all ranks in which B has "too many"
   if currpos[1][i] > B[i] then 
    currpos[1][i] := currpos[1][i] - 1; #B places card on pile
    currpos := switch(currpos);
    print(currpos);
    numberofturns:= numberofturns + 1;
    currpos[1][i] := currpos[1][i] + 1; #A takes card
    currpos := switch(currpos);
    print(currpos);
    numberofturns:= numberofturns + 1;
   else
    i := i + 1;
   fi;
  od;

  if sum(currpos[1][q], q = 1..k) >= 3 then  #if Player B now has at least 3 cards...
     
   for r from 1 to k while isAsturn = 0 do 
    if {seq(evalb(currpos[1][j] = 0),j = 1..r-1)} subset {true} and currpos[1][r] > 0 then 
     if currpos[1][r] >= 2 then #if Player B has at least two "lowest" ranking cards
      print("Switching turns:");
      currpos[1][r] := currpos[1][r] - 1; #1. B discards 'X'
      currpos := switch(currpos);
      print(currpos);
      numberofturns:= numberofturns + 1;
      currpos[1][r] := currpos[1][r] + 1; #2. A takes 'X'
      currpos := switch(currpos);
      print(currpos);
      numberofturns:= numberofturns + 1;
      currpos[1][r] := currpos[1][r] - 1; #3. B discards 'Y'
      currpos := switch(currpos);
      print(currpos);
      numberofturns:= numberofturns + 1;
      currpos[1][r] := currpos[1][r] - 1; #4. A discards 'x'
      currpos := switch(currpos);
      print(currpos);
      numberofturns:= numberofturns + 1;
      currpos[1][r] := currpos[1][r] + 2; #5. B takes 'X' and 'Y'
      currpos := switch(currpos);
      print(currpos);
      numberofturns:= numberofturns + 1;
      print("Turns switched.");
      isAsturn := 1;
     else 
      lowest := r;
     fi;
    elif lowest > 0 and currpos[1][r] > 0 then  #if Player B has a single lowest ranking card and r is the next highest one
      currpos[1][r] := currpos[1][r] - 1; #1. B discards 'X'
      currpos := switch(currpos);
      print("Switching turns:");
      print(currpos);
      numberofturns:= numberofturns + 1;
      currpos[1][r] := currpos[1][r] + 1; #2. A takes 'X'
      currpos := switch(currpos);
      print(currpos);
      numberofturns:= numberofturns + 1;
      currpos[1][lowest] := currpos[1][lowest] - 1; #3. B discards 'Y'
      currpos := switch(currpos);
      print(currpos);
      numberofturns:= numberofturns + 1;
      currpos[1][r] := currpos[1][r] - 1; #4. A discards 'x'
      currpos := switch(currpos);
      print(currpos);
      numberofturns:= numberofturns + 1;
      currpos[1][r] := currpos[1][r] + 1;     #5. B takes 'X' and 'Y'
      currpos[1][lowest] := currpos[1][lowest] + 1;
      currpos := switch(currpos);
      print(currpos);
      numberofturns:= numberofturns + 1;
      print("Turns switched.");
      isAsturn := 1;
    fi;
   od;

   #Now it's A's turn

   i := 1;
   while i <= k do   #keep going through all ranks in which A has "too many"
    if currpos[1][i] > A[i] then
     currpos[1][i] := currpos[1][i] - 1; #A places card on pile
     currpos := switch(currpos);
     print(currpos);
     numberofturns:= numberofturns + 1;
     currpos[1][i] := currpos[1][i] + 1; #B takes card
     currpos := switch(currpos);
     print(currpos);
     numberofturns:= numberofturns + 1;
    else
     i := i + 1;
    fi;
   od;

  else    #If Player B has fewer than 3 cards remaining
   print("Player B has fewer than 3 cards left.");
   for t from 1 to k do
    
    if currpos[1][t] > 0 and {seq(evalb(currpos[1][j]=0),j = 1..t-1)} subset {true} then #if i is smallest index for which B[i] > 0
      currpos[1][t] := currpos[1][t] - 1; #B places card on pile
      currpos := switch(currpos);
      print(currpos);
      numberofturns:= numberofturns + 1;
      currpos[1][t] := currpos[1][t] - 1; #A places same card
      currpos := switch(currpos);
      print(currpos);
      numberofturns:= numberofturns + 1;
      currpos[1][t] := currpos[1][t] + 2; #B takes
      currpos := switch(currpos);
      print(currpos);
      numberofturns:= numberofturns + 1;
      if B[t] = 4 then
       currpos[1][t] := currpos[1][t] - 1; #A places same card
       currpos := switch(currpos);
       print(currpos);
       numberofturns:= numberofturns + 1;
       currpos[1][t] := currpos[1][t] + 1; #B takes
       currpos := switch(currpos);
       print(currpos);
       numberofturns:= numberofturns + 1;
      fi;
    fi;
   od;
  fi;
 print("Number of turns taken so far:", numberofturns);

elif n > 0 then
 if (type(n, odd) or (n/2 + b < 3 or a < 1)) and (type(n,even) or (n-1)/2 + a < 3) then return("ERROR, violates conditions of corollary.") fi;
 currpos:= [A,B];
 highest := k;
 temppile:= [seq(0,i=1..k)]; #consists of the cards x1, x2, ..., xn from the corollary.
 while highest > 0 do
  if (highest = 1 and 3 - currpos[1][highest] - currpos[2][highest] > 0) or (highest > 1 and 4 - currpos[1][highest] - currpos[2][highest] > 0) then
   currpos := switch(currpos);
   currpos[1][highest] := currpos[1][highest] + 1;
   temppile[highest] := temppile[highest] + 1;
   ###print(temppile);
  else
   highest := highest - 1;
  fi;
 od;
 printpath(op(currpos));

 lowest := 1;

 while lowest <= k do
  if temppile[lowest] > 0 then
   temppile[lowest] := temppile[lowest] - 1;
   currpos[1][lowest] := currpos[1][lowest] - 1;
   currpos := switch(currpos);
   print(currpos);
  else
   lowest := lowest + 1;
  fi;
 od;
return;

fi;
  
end proc;


#attainable(PosSet): input a set of positions, outputs those which are attainable from a start position via the claim and corollary
attainable := proc(PosSet)
 local k, pos, a, b, n, i, j, V;
 k := nops(PosSet[1][1]);
 V := [];
 for pos in PosSet do
  a := sum(pos[1][i],i=1..k);
  b := sum(pos[2][j],j = 1..k);
  n := 4*k - 1 - a - b;
  if (type(n, odd) or (n/2 + b < 3 or a < 1)) and (type(n,even) or (n-1)/2 + a < 3) then
  else 
   V:= [op(V),pos];
  fi;
od;
return(V);
end proc;

#drawcomponent
drawcomponent:= proc(W,C)
 local V, i, G, H, legal, move, currpos, nextmove, cycle, j;
 G := [op({seq(W[i], i in C )} minus {-1, 1})];
 H:= [seq(op(i),i in G)]; #lists all drawn positions
 

 currpos:=H[1];
 V:=[[]];
 j:=1;
 while nops({op(H)} minus {seq(op(cycle), cycle in V)}) > 0 and j < 2 do  #while the union of the cycles do not form all of the drawn positions
 #print({seq(op(cycle), cycle in V)});
 
 while not member(currpos, V[j]) do
  
  V[j]:=[op(V[j]),currpos];
  legal:=legalmoves(op(currpos));
  for move in legal do

   if member(move, H) then  #if unassigned
    nextmove:=move;
    #print(nextmove);
   elif W[move] = -1 then   #make sure no child is assigned -1
    return("uh oh");
   end if;
  
  od;

  currpos:=nextmove;
    print(currpos);
 od;
 j:= j+1;
 
 #print( seq(op(cycle), cycle in V));
 #print(nops({op(H)} minus {seq(op(cycle), cycle in V)}));
 V:= [op(V),[]];
 #print(V);
 currpos:=({op(H)} minus {seq(op(cycle), cycle in V)})[1];

 od;
 return V;
end proc;

#legalmoves([2,2],[1,1]) should return [[[1,1],[1,0]], [[1,1],[1,2]]].  returns ordered n-tuple of resulting legal positions.
#Switches players A and B.

legalmoves := proc (A, B, rules::integer:=0) 
local S, i, j, maxcantake, l, newP, m, newA, k, highest, P, T, zero; 
k := nops(A);
zero := [seq(0,i=1..k)];
P := [seq(4, i = 1 .. k)]-A-B;
highest := 1;

if nops(B) <> k or P[1] < 1 or B = zero or not (member(false, {seq(evalb(A[i] = 0), i = 1 .. k)})) then  #error catching for input
 return []; 
else 
 T:=[];
 S := [];
 maxcantake := sum(P[i], i = 1 .. k)-1;



 if 0 < maxcantake then 	#step 2: check for 'take moves'
  newP := P;
  newA := A;
  for j to min(maxcantake, 3) do 
   l := k;
   while newP[l] = 0 and 1 <= l do 
    l := l-1;              #find the highest rank on the pile
   end do;
   if j = 1 then 
    highest := l;          #highest is the highest rank (index) on the original pile
   end if;
   newA[l] := newA[l]+1;    #take a card
   newP[l] := newP[l]-1;
  end do;
  S := [op(S), [B, newA]];
 end if;

for j from highest to k do	#step 1: check for 'placing' moves
 newA := A;
 newP := P;

 if j=1 and A[1] = 3 then	#step 1a: place 3 cards if they are all 9's
   newA[j] := 0;
   newP[j] := 4;
   T := [op(T), [B, newA]];
 elif A[j] <> 0 then	

  if A[j] = 4 then		#step 1b: place 4 cards
   newA[j] := 0;

   T := [op(T), [B, newA]];
   newA[j] :=3;
   T:= [op(T),[B,newA]];
  else
   newA[j] := newA[j]-1;         #step 1c: place one card
   T := [op(T), [B, newA]];
  end if;

 end if;

end do;

end if; 	 #end error catching
return [op(T),op(S)];   #put the placing moves at the beginning
end proc;


#uses old ordering of moves
oldlegalmoves := proc (A, B) 
local S, i, j, maxcantake, l, newP, m, newA, k, highest, P; 
k := nops(A);
P := [seq(4, i = 1 .. k)]-A-B;
highest := 1;

if nops(B) <> k or P[1] < 1 or not (member(false, {seq(evalb(A[i] = 0), i = 1 .. k)})) then  #error catching for input
 return []; 
else 
 S := [];
 maxcantake := sum(P[i], i = 1 .. k)-1;
 if 0 < maxcantake then 	#step 1: check for 'take moves'
  newP := P;
  newA := A;
  for j to min(maxcantake, 3) do 
   l := k;
   while newP[l] = 0 and 1 <= l do 
    l := l-1;
   end do;
   if j = 1 then 
    highest := l;
   end if;
   newA[l] := newA[l]+1;
   newP[l] := newP[l]-1;
  end do;
  S := [op(S), [B, newA]];
 end if;

for j from highest to k do	#step 2: check for 'placing' moves
 newA := A;
 newP := P;
 if A[j] <> 0 then		#step 2a: place one card
  newA[j] := newA[j]-1;
  newP[j] := newP[j]+1;
  S := [op(S), [B, newA]];
  if A[j] = 4 then		#step 2b: place 4 cards
   newA[j] := 0;
   newP[j] := 4;
   S := [op(S), [B, newA]]
  elif j = 1 and A[j] = 3 then	#step 2c: place 3 cards if they are all 9's
   newA[j] := 0;
   newP[j] := 4;
   S := [op(S), [B, newA]]
  end if;
 end if;
end do;
end if; 	 #end error catching
return S;
end proc;



#allpos(k): allpos(1) should return [[0],[0]],[[0],[1]],[[0],[2]],[[0],[3]],[[1],[0]],[[1],[1]],[[1],[2]]
allpos:=proc(k)  
local i,j,S,m,n,l, pos;  
S:=[];  
if k=1 then   
 S:=[seq(seq([[i],[j]],j=0..3-i),i=0..3)];    
  return(S);  
else  
 pos:=allpos(k-1);  
 for l from 1 to nops(pos) do    
  for m from 0 to 4 do    
   for n from 0 to 4-m do      
    S:=[op(S),[[op(pos[l][1]),m],[op(pos[l][2]),n]]];  
   od;  
  od;   
 od;  
fi;  
end proc;

#winorlose(k, moves): winorlose(2, 10) returns a table associating each position to either -1 or 1
winorlose := proc (k, moves) 
local Pos, n, i, j, wintable, l, m, newtable, p, legal; 
Pos := allpos(k); 
n := nops(Pos); 
wintable := table(); 
for i to moves do
 newtable := table(); 
 for j to n do 
  if not assigned(wintable[Pos[j]]) then
   if Pos[j][2] = [seq(0, l = 1 .. k)] then  #was Pos[j][1] = [0,0,...0] then label it a 1
    newtable[Pos[j]] := -1;
   else 
    legal := legalmoves(op(Pos[j])); 
    if {seq(wintable[legal[p]], p = 1 .. nops(legal))} = {1} then #if all children are W   
     newtable[Pos[j]] := -1; #call it L
    elif member(-1, {seq(wintable[legal[p]], p = 1 .. nops(legal))}) then #if a child is L     
     newtable[Pos[j]] := 1;  #call it W
    end if; 
   end if; 
  end if;
 end do; 
wintable := table([op(op(op(wintable))), op(op(op(newtable)))]);
end do; 
return wintable;
end proc;

#winlosemargin(k,intable,start,finish):
winlosemargin := proc (k, intable, start, finish)
local Pos, n, i, j, l, m, newtable, p, legal, moves, wintable;
wintable := intable;
moves := finish-start+1;
Pos := allpos(k);
n := nops(Pos); 
for i to moves do
 newtable := table();
 for j to n do
  if not assigned(wintable[j]) then 
   if Pos[j][1] = [seq(0, l = 1 .. k)] then 
    newtable[Pos[j]] := 1;
   else 
    legal := legalmoves(op(Pos[j])); 
    if {seq(wintable[legal[p]], p = 1 .. nops(legal))} = {1} then 
     newtable[Pos[j]] := -1;
    elif member(-1, {seq(wintable[legal[p]], p = 1 .. nops(legal))}) then 
     newtable[Pos[j]] := 1;
    end if;
   end if;
  end if; 
 end do; 
 wintable := table([op(op(op(wintable))), op(op(op(newtable)))]);
end do;
return wintable;
end proc;


#startpos(k): returns all starting positions for k-Pan.  startpos(1) returns [[2],[1]].  first coordinate goes through all possible ways to  have 2 k cards. 2 nd coordinate is determined.  
startpos:=proc(k) #is there a more efficient way?  `  
local i, S, n, Pos,j;   
Pos:=allpos(k);   
n:=nops(Pos);   
S:=[];  
for i from 1 to n do   
 if sum(Pos[i][1][j],j=1..k)=2*k and sum(Pos[i][2][j],j=1..k)=2*k-1 then    
  S:=[op(S),Pos[i]]   
 fi;  
od;
return(S);  
end proc;


#randpanmove(A,B): returns a random legal move from (A,B)
randpanmove:= proc(A,B)
 local legal, r;
 legal:=legalmoves(A,B);
 r:=rand(1..nops(legal));
 return(legal[r()]);
end proc;



#perfectpan(startA,startB): solves the game, spits out the winning move if it exists.
perfectpan:=proc(startA,startB)
local A,B,P, i, k, usermove, r, j, legal, depth, solvetable, q, p, legalcomp;
k:=nops(startA);
A:=startA;
B:=startB;
P:=[seq(4,j=1..k)]-A-B;
if k=2 then
 depth:=16;
elif k=3 then
 depth:=39;
fi;

solvetable:=winorlose(k,depth);

print("Pile:",P);
while B <> [seq(0,i=1..k)] and A <> [seq(0,i=1..k)] do
 legal:=legalmoves(A,B);
 if nops(legal) = 1 then
  legal:=[legal];
 fi;
 print("Choose a legal move:", seq(op(legal)[i][2],i=1..nops(legal) ));
 usermove:=readstat(terminal);
 legalcomp:=legalmoves(B,usermove);
 #TO DO: CHECK THAT THIS IS A LEGAL MOVE
 if solvetable[[B, usermove]][1]= 1 then
  print("It's looking good for me...");
  member(-1,[seq(solvetable[legalcomp[p]][1],p=1..nops(legalcomp))],'q'); #pick a move labelled with -1
  print([seq(solvetable[legalcomp[p]][1],p=1..nops(legalcomp))],member(-1,[seq(solvetable[legalcomp[p]][1],p=1..nops(legalcomp))]),q);
  A:=usermove;
  B:=legalcomp[q][2];

 elif solvetable[[B,usermove]][1]= -1 then
 print("Sigh...");
 member(solvetable[[B,usermove]][2],[seq(solvetable[legalcomp[p]][2],p=1..nops(legalcomp))],'q');
 
 print(solvetable[[B,usermove]][2], [seq(solvetable[legalcomp[p]][2],p=1..nops(legalcomp))], member(solvetable[[B,usermove]][2],[seq(solvetable[legalcomp[p]][2],p=1..nops(legalcomp))]),q);
 A:=usermove;
 B:=legalcomp[q][2];

 else 
 A:=usermove;
 #case for stalemates: check if any child is unassigned
 #B:=
 fi;

 P:=[seq(4,j=1..k)]-A-B;
 print("You chose: ",A,"Computer chose: ",B, "Pile:",P);
od;

if A= [seq(0,i=1..k)] then 
 print("You win!");
else
 print("Computer wins.");
fi;
end proc;


#winorlose3(k, moves): winorlose(2, 10) returns a table associating each position to either -1 or 1
#format: [[2,2],[1,1]]=[WL, nextpos, Length] = 
#[(1 or - 1 for W or L), (next best move), length of longest directed path]
winorlose3 := proc (k, moves) 
local Pos, n, i, j, wintable, l, m, newtable, p, legal, zed, pos2, x; 
zed:= [seq(0,l=1..k)];
Pos := allpos(k); 
n := nops(Pos); 
wintable := table(); 
for i to moves do
 newtable := table(); 
 for j to n do 
  if not assigned(wintable[Pos[j]]) then 
   if Pos[j][1] = zed then     #if player to move is the zero vector
    newtable[Pos[j]] := [1,Pos[j],1];
   else 
    legal := legalmoves(op(Pos[j]));
    pos2:=[seq(wintable[legal[p]][1], p = 1 .. nops(legal))];
    x:=0;
    if {op(pos2)} = {1} then #if all children are W   
     newtable[Pos[j]] := [-1,legal[1],i]; #call it L
    elif member(-1,pos2,'x') then #if a child is L  
     newtable[Pos[j]] := [1,legal[x],i];  #call it W
    end if; 
   end if; 
  end if;
 end do; 
wintable := table([op(op(op(wintable))), op(op(op(newtable)))]);
end do; 
return wintable;
end proc;




#evalscore(A,B): evalscore([2,2],[1,0]) will return the score associated with this position
#factor1: differences in amounts of low cards versus high cards
#factor3: points for having fewer cards
evalscore:=proc(A,B) #TO DO: also must keep in mind the number of groups of cards total.  smaller is better.  average per one extra k = 1.4
 local i,k,low, score, threesorfours, j, zed, cardratio, numgroupsA, numgroupsB, numcardsA, numcardsB;
 k := nops(A);
 zed:=[seq(0,j=1..k)];
 if k <> nops(B) then 
  return FAIL;
 end if;
 if A = zed then 
  return infinity;
 elif B = zed then
  return -infinity;
 else 
 numcardsA:=add(A[i], i= 1..k);
 numcardsB:=add(B[i], i= 1..k);
 low := (1-k)*(1/2);
 score := add((A[i]-B[i])*(low+i-1), i = 1 .. k);
#threesorfours := add(max(A[i]-3, 0)-max(B[i]-3, 0), i = 2 .. k)+max(A[1]-2, 0)-max(B[1]-2, 0); #Adds a point for each complete set you have, subtracts a point for each set your opponent has.
 numgroupsA:= numcardsA - add(4*max(A[i]-3, 0),i = 2..k) - 3*max(A[1]-2,0);
 numgroupsB:= numcardsB - add(4*max(B[i]-3, 0),i = 2..k) - 3*max(A[1]-2,0);
 cardratio:= (numgroupsB-numgroupsA)/2;
 return score + cardratio; #+threesorfours;
 end if;
end proc;


#minchild(A,B): returns a 2 element vector [loc, currmin] of the location and minimum score of each 'child' (legal move).  If A = 0, returns FAIL
minchild:= proc(A,B)
 local k,i,j,legal, currmin, loc, currscore;
 k:= nops(A);
 if k <> nops(B) then 
  return FAIL;
 end if;
 legal := legalmoves(A,B);
 currscore:=evalscore(op(legal[1]));
 currmin:=currscore;
 loc:=1;
 for j from 2 to nops(legal) do
  currscore:=evalscore(op(legal[j]));
  if currscore < currmin then 
   currmin:= currscore;
   loc:=j;
  end if;
 od;
 return [loc,currmin];
end proc;


#minchildpos #similar to minchild but returns the position as well which has the minimum score.
minchildpos:= proc(A,B)
 local k,i,j,legal, currmin, loc, currscore, zed, l;
 k:= nops(A);
 zed:=[seq(0,l=1..k)];
 if k <> nops(B) then 
  return FAIL;
 end if;
 legal := legalmoves(A,B);
 currscore:=evalscore(op(legal[1])); 
 currmin:=currscore;
 loc:=1;
 for j from 2 to nops(legal) do         ##can make this more efficient by stopping the search when you find -infinity
  currscore:=evalscore(op(legal[j]));
  if currscore < currmin then 
   currmin:= currscore;
   loc:=j;
  end if;
 od;
 return [loc,currmin,legal[loc]];
end proc;

#minmaxonestep: returns a 2 element vector [loc, currmax] of the location and value of the child whose children have the largest minimum score
minmaxonestep:= proc(A,B)
 local i, k, legal, currmax, loc, currscore, j;
 k:= nops(A);
 if k <> nops(B) then
  return FAIL;
 end if;
 legal := legalmoves(A,B);
 currscore:=minchild(op(legal[1]))[2];
 currmax:=currscore;
 loc:=1;
 for j from 2 to nops(legal) do
  currscore:=minchild(op(legal[j]))[2];
  if currscore > currmax then 
   currmax:= currscore;
   loc:=j;
  end if;
 od;
 return [loc,currmax];

end proc;

#minmaxonesteppos: similar to minmaxonestep, but returns the position as well
minmaxonesteppos:= proc(A,B)
 local i, k, legal, currmax, currscore, j, currchild, currpos, loc;
 k:= nops(A);
 if k <> nops(B) then
  print(A,B);
  return FAIL;
 end if;
 legal := legalmoves(A,B);
 currchild:=minchildpos(op(legal[1]));
 currscore:=currchild[2];
 currpos:=currchild[3];
 currmax:=currscore;
 loc:=1;
 for j from 2 to nops(legal) do              #for each CHILD
  currchild:=minchildpos(op(legal[j]));      #find the minimum grandchild
  currscore:=currchild[2];
  if currscore > currmax then                #take the maximum of these minima
   currmax:= currscore;
   currpos:=currchild[3];
   loc:=j;
  end if;
 od;
 return [loc,currmax, currpos];  #return recommended move number, score for that move, resulting worst case scenario position

end proc;



#returns a 3 element vector of the location and value and position associated with the child with the iterated minmax with depth
minmaxpan:= proc(A,B, depth)
 local i;
  if depth = 1 then
   return minmaxonesteppos(A,B);
  elif depth > 1 then
    return minmaxpan(op(minmaxonesteppos(A,B)[3]),depth-1);  #need to send 2 and 3
  end if;

end proc;


minmaxrecursive:= proc(A,B,score)
 local i, k, legal, currmax, currscore, j, currchild, currpos, loc;
 k:= nops(A);
 if k <> nops(B) then
  print(A,B);
  return FAIL;
 end if;
 legal := legalmoves(A,B);
 currchild:=minchildpos(op(legal[1]));    ###can make this more efficient by stopping the search when you find +infinity
 currscore:=currchild[2];
 currpos:=currchild[3];
 currmax:=currscore;
 loc:=1;
 for j from 2 to nops(legal) do              #for each CHILD
  currchild:=minchildpos(op(legal[j]));      #find the minimum grandchild
  currscore:=currchild[2];
  if currscore > currmax then                #take the maximum of these minima
   currmax:= currscore;
   currpos:=currchild[3];
   loc:=j;
  end if;
 od;
 return [loc,currmax, currpos];  #return recommended move number, score for that move, resulting worst case scenario position


end proc;



#switch(Pos): Switches the turn order.  Ex: switch([[2,2],[1,2]]) returns [[1,2],[2,2]]
switch:= proc(Pos)
 local newPos;
 newPos:= [Pos[2],Pos[1]];
 return newPos;
end proc;

#adapted from pseudocode from Wikipedia article on minimax algorithm
#TO DO: return principal variation
minmaxfinal:=proc(A,B, depth, maximizingPlayer)
 local bestvalue, legal, child, val, princvar, bestchild, multiplier;
 legal:=legalmoves(A,B);
 multiplier:= -1;
 if maximizingPlayer then multiplier:=1 end if;

 if depth <= 0 or nops(legal) = 0 then
#   return evalscore(A,B); 
  return [multiplier*evalscore(A,B),[]];
 end if;

 if maximizingPlayer then 
  multiplier = 1;    
  bestvalue:= -infinity;
  bestchild:= legal[1];
 for child in legal do
  val:=minmaxfinal(op(child), depth - 1, false)[1];
# print(val);
  if val > bestvalue then               #update current maximum of all children
    bestvalue:=val;
    bestchild:=child;
  end if;
  princvar:= [bestchild, op(princvar)];
 od;
 return [bestvalue, bestchild];

 else
  bestvalue:= infinity;
  bestchild:= legal[1];
  for child in legal do
   val:=minmaxfinal(op(child), depth - 1, true)[1];
   if val < bestvalue then
    bestvalue:=val;                #update current minimum of all children
    bestchild:=child;
   end if;
   princvar:= [bestchild, op(princvar)];
  od;
  return [bestvalue, bestchild];
 end if;

end proc;


minimaxpan:=proc(startA,startB, depth, advisor::integer:=0)
local zed, A,B,P, i, k, usermove, legal, legalcomp, legit, wintable, prevA, prevB, movehist, fakedepth, Bvector;
fakedepth:=depth;
k:=nops(startA);
zed:=[seq(0,i=1..k)];
A:=startA;
B:=startB;
P:=[seq(4,j=1..k)]-A-B;
if advisor>0 then
 print("Loading PerfectPan, please wait...");
 wintable:=winorlose3(k,advisor);
fi;

print("Pile:",P);
while B <> zed and A <> zed do
 legal:=legalmoves(A,B);
 legit:=seq(op(legal[i])[2],i=1..nops(legal));
 do
  print("Choose a legal move:",legit);
  if advisor>0 then
   print("PerfectPan advises",wintable[[A,B]]);
  fi;
  usermove:=readstat(terminal);
  if member( usermove, {legit} ) then 
   break;
  elif usermove= "undo" then
   A:=prevA;
   B:=prevB;
   break;
  else print("Sorry, I didn't catch that.");
  fi;
 od;
 if usermove <> "undo" then
  prevA:=A;
  prevB:=B;
  A:=usermove;
  do
   Bvector:=minmaxfinal(B,usermove,fakedepth,true);
   if Bvector[1] = -infinity and fakedepth> 2 then
    fakedepth:=fakedepth-2;    #used in order to decide when all moves lead to sure defeat, choose the longest path.
   else
    B:= Bvector[2][2];
    break;
   fi;
  od;
 end if;

 P:=[seq(4,j=1..k)]-A-B;
 print("You chose: ",A,"Computer chose: ",B, "\nPile:",P);
od;

if A= zed then 
 print("You win!");
else
 print("Computer wins.");
fi;
end proc;


#randstart(k): returns a random starting position from a k-Pan game
#does this by listing cards from 1 to 4k-1.  chooses 2k-1 distinct numbers (this is for the guy who just moved)
randstart:=proc(k)
 local i, numcards, j, allnums, currnums, randpos, choosenum, chosennums, thispos, chosenvector, B,A;
 allnums:=[seq(j,j=1..4*k-1)];
 currnums:=allnums;
 chosennums:=[];
 for i from 1 to 2*k-1 do
  randpos:=rand(1..nops(currnums));
  thispos:=randpos();
  chosennums:=[op(chosennums),currnums[thispos]];
  currnums:=remove(j->j=currnums[thispos],currnums);
 od;

 chosenvector:=(chosennums-(chosennums mod 4))/4;
 B:=[seq(0,j=1..k)];
 for i from 1 to nops(chosenvector) do
  B[chosenvector[i]+1]:=B[chosenvector[i]+1]+1;     #correction to get the index right.  bumps up the entry in the chosenvector[i]+1'th component
 od;
 A:=[3,seq(4,j=2..k)]-B;
 return(A,B);
end proc;


#alphabeta(A,B,depth,alpha,beta,maximizingPlayer): alpha beta pruning search for best moves
alphabeta:=proc(A,B,depth, alpha, beta, maximizingPlayer)
 local multiplier, legal, child, tempalpha, tempbeta, val, bestval, bestchild;
 tempalpha:=alpha;
 tempbeta:=beta;
 multiplier:= -1;
 if maximizingPlayer then multiplier:=1 end if;
 legal:=legalmoves(A,B);
 if depth <= 0 or nops(legal) = 0 then
   return [multiplier*evalscore(A,B),[]];
 end if;

  bestchild:=legal[1]; 

 if maximizingPlayer then
 
   bestval:= -infinity;
   for child in legal do
    val:= alphabeta(op(child), depth - 1, tempalpha, tempbeta, false)[1];
    if tempalpha < val then
    #tempalpha := max(tempalpha, alphabeta(op(child), depth - 1, tempalpha, tempbeta, false));
     bestchild:=child;
     bestval:=val;
     tempalpha:=val;
     if tempbeta <= tempalpha then
      break;    # beta cut-off
     end if;
    end if;
   od;
   return [tempalpha,bestchild];
 else
   for child in legal do
    bestval:=infinity;
    val:= alphabeta(op(child), depth - 1, tempalpha, tempbeta, true)[1];
    if tempbeta > val then
    bestchild:=child;
    bestval:=val;
    tempbeta:=val;
    #tempbeta := min(tempbeta, alphabeta(op(child), depth - 1, tempalpha, tempbeta, true));
    if tempbeta <= alpha then
     break;  # alpha cut-off
    end if;
    end if;
   od;
   return [tempbeta,bestchild];
 end if;
end proc;

#alphabetapan uses alpha beta pruning
alphabetapan:=proc(startA,startB, depth, advisor::integer:=0)
local zed, A,B,P, i, k, usermove, legal, legalcomp, legit, wintable, prevA, prevB, movehist, fakedepth, Bvector, advisormove;

k:=nops(startA);
zed:=[seq(0,i=1..k)];
A:=startA;
B:=startB;
P:=[seq(4,j=1..k)]-A-B;
if advisor=2 then
 print("Loading PerfectPan, please wait...");
 wintable:=winorlose3(k,39);
fi;

print("Pile:",P);
while B <> zed and A <> zed do
 legal:=legalmoves(A,B);
 legit:=seq(op(legal[i])[2],i=1..nops(legal));
 do
  print("Choose a legal move:",legit,"undo");
  if advisor=2 then
   print("PerfectPan advises",wintable[[A,B]]);
  elif advisor=1 then
   advisormove:=alphabeta(A,B,depth,-infinity,infinity,true);
   print("alphabetapan advises",advisormove[2][2],"which has a score of",advisormove[1]);
  fi;
  usermove:=readstat(terminal);
  if member( usermove, {legit} ) then 
   break;
  elif usermove= "undo" then
   A:=prevA;
   B:=prevB;
   break;
  else print("Sorry, I didn't catch that.");
  fi;
 od;
 if usermove <> "undo" then
  prevA:=A;
  prevB:=B;
  A:=usermove;
  fakedepth:=depth;
  do
   Bvector:=alphabeta(B,usermove,fakedepth,-infinity,infinity,true);
   if Bvector[1] = -infinity and fakedepth> 2 then
    fakedepth:=fakedepth-2;    #used in order to decide when all moves lead to sure defeat, decrease the depth until a best move is found.
   else
    B:= Bvector[2][2];
    break;
   fi;
  od;
 end if;

 P:=[seq(4,j=1..k)]-A-B;
 print("You chose: ",A,"Computer chose: ",B, "\nPile:",P);
od;

if A= zed then 
 print("You win!");
else
 print("Computer wins.");
fi;
end proc;


#randpos(k): returns a random position from a k-Pan game.
#This procedure does this by choosing a random solution to a + b + p = n for each rank, where a, b are nonnegative integers, 
#(n = 3 and p is positive) if the rank is 9, and (n = 4 and p is nonnegative) otherwise.  
randpos:=proc(k)
 local i, j, containersetnine, containersetother, choose1, choose2, B, A, P; 
 

 choose1:=rand(1..5);
 choose1:=choose1();
 containersetnine:=[op({seq(i, i = 1..5)} minus {choose1})];    #to solve, x1 + ...+ xm = n for xi nonnegative, line up n+m-1 containers, choose m-1 of them to be "+" signs
 choose2:=rand(1..4);
 choose2:=choose2();
 choose2:=containersetnine[choose2];  #choose a number distinct from the first
 A:= [min(choose1,choose2)-1];      #These three lines translate the placement of the "+" signs to a solution in nonnegative integers
 B:= [max(choose1,choose2)-A[1]-2];  
 #P:= [4-A[1]-B[1]];

 for j from 2 to k do
  choose1:=rand(1..6);
  choose1:=choose1();
  containersetother:=[op({seq(i, i = 1..6)} minus {choose1})];    #do the same thing for the ranks other than 9
  choose2:=rand(1..5);
  choose2:=choose2();
  choose2:=containersetother[choose2];
  A:= [op(A),min(choose1,choose2)-1];      
  B:= [op(B),max(choose1,choose2)-A[j]-2];  
  #P:= [op(P),4-A[j]-B[j]];
 od;

 return(A,B);
end proc;

#picks a random sample of n k-pan hands, returns an array of these hands
randsamp:=proc(k,n)
 local S;
 S:= {};
 while nops(S) < n do;
  S:= S union {[randpos(k)]};
 od;
 return [op(S)];

end proc;