#######################################################################
## Research: Save this file as LRC.txt. To use it, stay in the        #
## same directory, get into Maple (by typing: maple <Enter> )         #
## and then type:  read "LRC.txt"; <Enter>                            #
## Then follow the instructions given there                           #
##                                                                    #
## Written by Bryan Ek, Rutgers University,                           #
##  bryan.t.ek@math.rutgers.edu.                                      # 
#######################################################################

print(`This is LRC`):
print(`by Bryan Ek, advisee to Dr. Zeilberger`):
print(``):
print(`The latest version (5/1/2017) of this package is available from`):
print(`http://www.math.rutgers.edu/~bte14/Code/LeftRightCenter/LRC.txt`):
print(``):
print(`Left Right Center is a completely random dice game involving moving chips amongst players based on dice rolls.`):
print(`As there is no strategy to affect the game, and each state can lead to only finitely many other states,`):
print(`we can find the expected length of a game by enumerating all expressions and then using a computer to solve the (usually very large) linear system of equations.`):
print(``):
print(`For a list of functions this package provides type Help().`):
print(``):

with(combinat):

Help:=proc():
	if args=`ExpectedLengthEquation` then
		print(`ExpectedLengthEquation(L,E,maxDice): Given a starting position, creates an equation for how long the game is expected to last.`):
		print(`L is a list of non-negative integers. E is a symbol. maxDice is a positive integer.`):
		print(`Try ExpectedLengthEquation([3,2],E,3).`):
	elif args=`ExpectedLength` then
		print(`ExpectedLength(numChips,numPlayers,E,maxDice): Given a starting position, creates all equations and solves for how long the game is expected to last for any distribution of numChips.`):
		print(`numChips is a non-negative integer. numPlayers is a positive integer. E is a symbol. maxDice is a positive integer.`):
	elif args=`DistributeChips` then
		print(`DistributeChips(numChips,numPlayers): Outputs all possible distributions of numChips amongst numPlayers.`):
	else
		print(`Current useable functions:`):
		print(`ExpectedLengthEquation, ExpectedLength, DistributeChips`):
		print(`For help with a specific function, type Help(function).`):
	fi:
end:

######################################################################################################
#ExpectedLengthEquation(L,E,maxDice): Given a starting position, creates an equation for how long the game is expected to last.
#L is a list of non-negative integers. E is a symbol. maxDice is a positive integer.
ExpectedLengthEquation:=proc(L,E,maxDice) local n,numDiceToRoll:
	n:=nops(L):
	if n<2 or convert(L,`+`)=max(L) then#Either there are not enough people to play, or only 1 person has chips remaining.
		RETURN(E[op(L)]=0):
	fi:
	if L[1]=0 then
		E[op(L)]=E[op(2..n,L),L[1]]:#A roll never takes place. This could be edited to count null turns.
	elif n=2 then #Left and right rolls are equivalent.
		numDiceToRoll:=min(maxDice,L[1]):
		E[op(L)]=add(add((1+E[L[2]+a,L[1]-a-c])*multinomial(numDiceToRoll,a,c,numDiceToRoll-a-c)/(3^(a+c)*2^(numDiceToRoll-a)),c=0..numDiceToRoll-a),a=0..numDiceToRoll):
	else #n>2.
		numDiceToRoll:=min(maxDice,L[1]):
		E[op(L)]=add(add(add((1+E[L[2]+a,op(3..n-1,L),L[n]+b,L[1]-a-b-c])*multinomial(numDiceToRoll,a,b,c,numDiceToRoll-a-b-c)/(3^(a+b+c)*2^numDiceToRoll),c=0..numDiceToRoll-a-b),b=0..numDiceToRoll-a),a=0..numDiceToRoll):
	fi:
end:

###Then run through, put stuff into set that has already found equations, keep going until toFind set is empty.
###Can ignore anything that has fewer chips since this could have already been solved for.

######################################################################################################
#ExpectedLength(numChips,numPlayers,E,maxDice): Given a starting position, creates all equations and solves for how long the game is expected to last for any distribution of numChips.
#numChips is a non-negative integer. numPlayers is a positive integer. E is a symbol. maxDice is a positive integer.
ExpectedLength:=proc(numChips,numPlayers,E,maxDice) option remember:
	if numChips=1 or numPlayers=1 then
		return(ExpectedLength1(numChips,numPlayers,E,maxDice)):
	fi:
	subs({seq(op(ExpectedLength(numChips-i,numPlayers,E,maxDice)),i=1..min(numChips-1,maxDice))},ExpectedLength1(numChips,numPlayers,E,maxDice)):
end:
######################################################################################################
#ExpectedLength1(numChips,numPlayers,E,maxDice): Given a starting position, creates all equations and solves for how long the game is expected to last for any distribution of numChips.
#numChips is a non-negative integer. numPlayers is a positive integer. E is a symbol. maxDice is a positive integer.
ExpectedLength1:=proc(numChips,numPlayers,E,maxDice) local distributions,equations:
	distributions:=DistributeChips(numChips,numPlayers):
	equations:={seq(ExpectedLengthEquation(L,E,maxDice),L in distributions)}:
	solve(equations,{seq(E[op(L)],L in distributions)}):
end:


######################################################################################################
#DistributeChips(numChips,numPlayers): Outputs all possible distributions of numChips amongst numPlayers.
DistributeChips:=proc(numChips,numPlayers) option remember:
	if numChips=0 then
		return({[0$numPlayers]}):
	elif numPlayers=1 then
		return({[numChips]}):
	fi:
	{seq(seq([i,op(L)],L in DistributeChips(numChips-i,numPlayers-1)),i=0..numChips)}:
end: