######################################################################
##DiceGoals.txt: Save this file as DiceGoals.txt.                    #
#To use it, stay in the                                              #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read `DiceGoals.txt` : <Enter                      #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Lucy Martinez with help from Doron Zeilberger           #
######################################################################
 
#Version: February 2026
with(Statistics):


print(`First Written: November 2025: tested for Maple 20  `):
print():
print(`This is DiceGoals.txt, a Maple package, accompanying the article:`):
print(`How Many Dice Rolls Does It Take to Achieve Your Prescribed Goal?`):
print(`by Christoph Koutschan and Lucy Martinez`):


print(`This Maple package, as well as the article, are `):
print(`available from:`):
print(` WEBSITE HERE`):
print(``):
print(`Please report bugs to lucy.martinez@rutgers.edu `):

print(``):
print(`-----------------------------`):

print(``):
print(`For a list of the procedures type Help(); `):
print(`for help with a specific procedure, type Help(procedure_name)`):
print(``):
print(`For a list of the supporting functions type: Help1();`):
print(`-----------------------------`):
print(``):
print(``):


print(`For a list of the recurrence package provided by Zeilberger type HelpRec(); `):
print(`for help with a specific procedure, type HelpRec(procedure_name)`):
print(`-----------------------------`):
print(``):
print(``):


print(`For a list of the faster procedures and recurrences type HelpFast(); `):
print(`for help with a specific procedure, type Help(procedure_name)`):
print(`-----------------------------`):
print(``):
print(``):


print(`For a list of the simulation procedures type HelpSimu(); `):
print(`for help with a specific procedure, type Help(procedure_name)`):
print(`-----------------------------`):
print(``):
print(``):


#print(`For a list of the Story functions type:  HelpS();`):
#print(`for help with a specific procedure, type Help(procedure_name); `):
#print(`-------------------------------`):
#print(``):
#print(``):




HelpSimu:=proc()
if args=NULL then
print(`The simulation procedures are: easySimuG, generalSimuG, ManyeasySimuG, ManygeneralSimuG, ManySimuG, ManyvgSimuG,`):
print(`Simu, SimuG, vgSimuG`):
else
Help(args):
fi:
end:

HelpS:=proc()
if args=NULL then

print(`The Story procedures are: `):

else
 Help(args):
fi:
end:

Help1:=proc()
if args=NULL then
print(`The supporting functions are: CheckEasyGenderRatioConjecture, CheckGenderRatioConjecture, CheckGenderRatioConjectureVG,`):
print(`CheckP, Chop, comp, easyChop, easycomp, easyGFB, easyGFBd, easyGFBt, generalcomp, generalChop,`):
print(`RandomLinear, RandomS, Seqs, Seqs1, vgChop, vgcomp`):
else
Help(args):
fi:
end:

HelpFast:=proc()
if args=NULL then
print(`The faster procedures are: Ck, ExpF, ExpG, ExpFgen, F, Fgen, G`):
print(`OpeAve2, Ope3same, Ope4same, Ope5same, Ope6same`):
print(`Seq2same,Seq3same,Seq4same,Seq5same,Seq6same`):
else
Help(args):
fi:
end:



Help:=proc()

if args=NULL then
print(`The main procedures are: GFB, GFBC, GFBCt, GFBCgeneral, GFBgeneral, GFBCgeneralT,`):
print(`GFBCvg, GFBvg, GFBvgT,`):
print(`Stat, StatEasyFS, StatFS`):
print(`vgExpAND, vgExpOR, vgPGAND, vgPGOR`):
#Trend

elif nargs=1 and args[1]=CheckEasyGenderRatioConjecture then 
print(`CheckEasyGenderRatioConjecture(G,P): Inputs a list G corresponding to the number of each genders, and`):
print(`a list of probabilities P and outputs the ratio of the expected babies from each gender`):
print(`and their corresponding probability such that there is only one of the genders reached as a goal. Try: `):
print(`CheckEasyGenderRatioConjecture([5,6,8],[1/6,1/3,1/2]); `):


elif nargs=1 and args[1]=CheckGenderRatioConjectureVG then 
print(`CheckGenderRatioConjectureVG(S,P,K): Inputs a set S where a member Si of S1`):
print(`is a set of functions in x[1],...,x[k] and k is the size of P`):
print(`a list of probabilities P and an integer K and outputs the ratio of the expected babies from each gender`):
print(`and their corresponding probability for the general case in which exactly one set in S is non-negative. Try:`):
print(`CheckGenderRatioConjectureVG({{x[1]-5,x[2]-6,x[3]-8}},[1/6,1/3,1/2],100); `):


elif nargs=1 and args[1]=CheckGenderRatioConjecture then 
print(`CheckGenderRatioConjecture(G,P,K): Inputs a list G corresponding to the number of each genders,`):
print(`a list of probabilities P and an integer K and outputs the ratio of the expected babies from each gender`):
print(`and their corresponding probability for the general case in which every goal is met. Try: `):
print(`CheckGenderRatioConjecture([5,6,8],[1/6,1/3,1/2],100); `):


elif nargs=1 and args[1]=CheckP then 
print(`CheckP(P): Given a probability table outputs true if the entries`):
print(`add up to 1. Try: `):
print(`CheckP([1/2,1/4,1/4]); `):

 
elif nargs=1 and args[1]=Chop then 
print(`Chop(F,x,G): given a polynomial F in x[1],...,x[k] where k is the number of faces on the die,`):
print(`and a list G with the number of (desired) occurrences of each side of a die, outputs the part with exponents`):
print(`that satisfy reaching the goal of having each of the sides of the die in the list G. Try: `):
print(`Chop(1/5*x[1]^3*x[2]^8 + 1/5*x[1]^3*x[2]^7 ,x,[3,8]); `):

 
elif nargs=1 and args[1]=Ck then 
print(`Ck(k,K): inputs integers k and K where 3<=k<=6 and outputs the constant Ck where`):
print(`Ck is asymptotically equal to (k*N-ak(K))/(k*sqrt(K)) as K goes to infinity`):
print(`where ak(K) is the expected number of die rolls of a k-sided die such that`):
print(`your goal is to see K ocurrences of one side and each side has prob. 1/k. Try:`):
print(`Ck(3,1000); `):


elif nargs=1 and args[1]=comp then 
print(`comp(L1,L2): Inputs two lists L1 and L2 and compares their values`):
print(`outputs true if the elements in L2 are greater or equal to`):
print(`the elements in L1 (pointwise). Try: `):
print(`comp([2,2],[4,3]); `):


elif nargs=1 and args[1]=easyChop then 
print(`easyChop(F,x,G): given a polynomial F in x[1],...,x[k] where k is the number of faces on the die,`):
print(`and a list G with the number of (desired) occurrences of each side of a die, outputs the part with exponents`):
print(`that satisfy reaching the goal of reaching all the occurrences of one of the sides of the die from the list G. Try: `):
print(`easyChop(1/5*x[1]^3*x[2]^8 + 1/5*x[1]^3*x[2]^7 ,x,[9,8]); `):


elif nargs=1 and args[1]=easycomp then 
print(`easycomp(L1,L2): Inputs two lists L1 and L2 and compares their values`):
print(`outputs true if one of the elements in L2 are greater or equal to`):
print(`the corresponding entry in L1 (pointwise). Try: `):
print(`easycomp([2,4],[2,1]); `):


elif nargs=1 and args[1]=easyGFB then 
print(`easyGFB(G,P,x): inputs a list G corresponding to the number of (desired) occurrences of each side of a die,`):
print(`a list of probabilities P (for each side of the die), a variable x,`):
print(`outputs the prob. generating function in x[1],...,x[k] where k is the number of faces on the die`):
print(`whose coefficient of x[1]^a[1]*...*x[k]^a[k] is the prob. that once you have reached your goal the prob.`):
print(`either a[1]=g[1] (if you reached G[1]) or a[2]=g[2] (if you reached G[2]), etc. Try:`):
print(`easyGFB([5,3,2],[1/2,1/8,3/8],x); `):


elif nargs=1 and args[1]=easyGFBd then 
print(`easyGFBd(G,P,x): Same as easyGFB(G,P,x) (q.v), but faster, since it is done by direct summation. Try:`):
print(`easyGFBd([5,3,2],[1/2,1/8,3/8],x); `):

elif nargs=1 and args[1]=easyGFBt then 
print(`easyGFBt(G,P,t): It computes easyGFB(G,P,x) and substitutes each variable`):
print(`x[1],x[2],...,x[k] (k is the size of G) to one variable t indicating the total number of dice rolls. Try:`):
print(`easyGFBt([5,3,2],[1/2,1/8,3/8],t); `):


elif nargs=1 and args[1]=easySimuG then 
print(`easySimuG(L,P): Inputs a list of positive integers L, and`):
print(`a list of non-neg. rational numbers P that adds up 1, of the same length as L, and`):
print(`outputs a list X where there exists one entry in X`):
print(`that is greater or equal to the corresponding value in L (pointwise)`):
print(`where X was first a list of zeros and with prob. P[i],`):
print(`X[i] was increased by 1 each time. Try:`):
print(`easySimuG([3,4,2],[1/3,1/3,1/3]); `):

 
elif nargs=1 and args[1]=ExpF then 
print(`ExpF(n,p1,p2,p3): outputs the first k terms of the expected duration of the game`):
print(`in which with prob. p1 the die lands on side 1, with prob. p2 the die lands on side 2,`):
print(`with prob. p3 the die lands on side 3, and`):
print(`at each die roll, you keep track of the number of appearances of each side,`):
print(`the game stops until either the number of appearances of side 1 equals k,`):
print(`OR the number of appearances of side 2 equals k,`):
print(`OR the number of appearances of side 3 equals k, where 1<=k<=n. Try:`):
print(`ExpF(10,1/3,1/3,1/3); `):

 
elif nargs=1 and args[1]=ExpFgen then 
print(`ExpFgen(n,P): outputs the first n terms of the expected duration of the game`):
print(`in which with probability P[i] the die lands on side i for 1<=i<=k for each k<=n`):
print(`and at each die roll, you keep track of the number of appearances of each side,`):
print(`the game stops until either the number of appearances of side 1 equals k,`):
print(`OR the number of appearances of side 2 equals k,`):
print(`OR the number of appearances of side 3 equals k,`):
print(`OR the number of appearances of side 4 equals k, where 1<=k<=n. Try:`):
print(`ExpFgen(10,[1/3,1/3,1/3]); `):


elif nargs=1 and args[1]=ExpG then 
print(`ExpG(n,p1,p2,p3,p4): outputs the first n terms of the expected duration of the game`):
print(`in which with prob. p1 the die lands on side 1, with prob. p2 the die lands on side 2,`):
print(`with prob. p3 the die lands on side 3, with prob. p4 the die lands on side 4, and`):
print(`at each die roll, you keep track of the number of appearances of each side,`):
print(`the game stops until either the number of appearances of side 1 equals k,`):
print(`OR the number of appearances of side 2 equals k,`):
print(`OR the number of appearances of side 3 equals k,`):
print(`OR the number of appearances of side 4 equals k, where 1<=k<=n. Try:`):
print(`ExpG(10,1/4,1/4,1/4,1/4);`):


elif nargs=1 and args[1]=F then 
print(`F(m1,m2,m3,p1,p2,p2,t): Inputs non-negative integers m1, m2, and m3,`):
print(`probabilities p1,p2,p3, and a variable t, outputs`):
print(`the probability generating function of the duration of the game in which`):
print(`with prob. p1 the die lands on side 1, with prob. p2 the die lands on side 2,`):
print(`with prob. p3 the die lands on side 3 and`):
print(`at each die roll, you keep track of the number of appearances of each side,`):
print(`the game stops until either the number of appearances of side 1 equals m1,`):
print(`OR the number of appearances of side 2 equals m2, OR`):
print(`the number of appearances of side 3 equals m3. Try:`):
print(`F(5,5,5,1/3,1/3,1/3,t); `):


elif nargs=1 and args[1]=Fgen then 
print(`Fgen(M,P,t): inputs a list M and a list of probabilities P, each of length k,`):
print(`and a variable t, outputs the probability generating function of the duration of the game in which`):
print(`with probability P[i] the die lands on side i for 1<=i<=k and the game stops until`):
print(`the number of appearances of side 1 equals M[1]`):
print(`OR the number of appearances of side 2 equals M[2],`):
print(`OR the number of appearances of side 3 equals M[3],...,`):
print(`OR the number of appearances of side k equals M[k]. Try:`):
print(`Fgen([3,3,3],[1/3,1/3,1/3],t);`):


elif nargs=1 and args[1]=G then 
print(`G(m1,m2,m3,m4,p1,p2,p2,p4,t): Inputs non-negative integers m1, m2, m3, and m4`):
print(`probabilities p1,p2,p3,p4, and a variable t, outputs the`):
print(`probability generating function of the duration of the game in which`):
print(`with prob. p1 the die lands on side 1, with prob. p2 the die lands on side 2,`):
print(`with prob. p3 the die lands on side 3, with prob. p4 the die lands on side 4,`):
print(`and at each die roll, you keep track of the number of appearances of each side,`):
print(`the game stops until either the number of appearances of side 1 equals m1,`):
print(`OR the number of appearances of side 2 equals m2,`):
print(`OR the number of appearances of side 3 equals m3,`):
print(`OR the number of appearances of side 4 equals m4. Try:`):
print(`G(5,5,5,5,1/4,1/4,1/4,1/4,t); `):


elif nargs=1 and args[1]=generalChop then 
print(`generalChop(F,x,G,g): given a polynomial F in x[1],...,x[k] where k is the number of faces on the die,`):
print(`and a list G with number of (desired) occurrences of each side of a die, outputs the part with exponents`):
print(`that satisfy reaching the goal of having exactly g out of the nops(G) desired ocurrences from the list G. Try:`):
print(`generalChop(1/5*x[1]^3*x[2]^8 + 1/5*x[1]^3*x[2]^7,x,[3,8],2); `):


elif nargs=1 and args[1]=generalcomp then 
print(`generalcomp(L1,L2,g): inputs two lists L1 and L2 (same length) and an integer g<=nops(L1)`):
print(`compares their values and outputs true if there are EXACTLY g elements in L1 and L2 such that`):
print(`those elements in L2 are greater or equal to the corresponding entry in L1 (pointwise). Try: `):
print(`generalcomp([1,2,3],[2,0,5],2); `):


elif nargs=1 and args[1]=generalSimuG then 
print(`generalSimuG(L,P,g): inputs a list of positive integers L,`):
print(`a list of non-neg. rational numbers P that adds up 1, of the same length as L, and`):
print(`an integer g<=nops(L) and outputs a list X where there exists EXACTLY g entries in X`):
print(`that are greater or equal to the corresponding value in L (pointwise)`):
print(`where X was first a list of zeros and with prob. P[i],`):
print(`X[i] was increased by 1 each time. Try: `):
print(`generalSimuG([3,2,2],[1/3,1/3,1/3],1); `):


elif nargs=1 and args[1]=GFB then 
print(`GFB(G,P,x, K): Inputs a list G corresponding to the number of each desired occurence`):
print(`of each side on a die with k:=nops(G) faces, a list of probabilities P (for each side of the die),`):
print(`a variable x, and an integer K, outputs the prob. generating function in x[1],...,x[k]`):
print(`whose coefficient of x[1]^a[1]*...*x[k]^a[k] is the prob. that once you have reached your goal the prob.`):
print(`a[i]=g[i] (when you reached G[i] last) for some i, and for all j not equal to i, a[j]>=g[j] (when you reach G[j] before reaching G[i]),`):
print(`assuming that you finished in <=K steps, in other words, you reach all your goals,`):
print(`and also outputs the probability of not reaching your goal. Try:`):
print(`GFB([5,3,2],[1/2,1/8,3/8],x,10); `):
 

elif nargs=1 and args[1]=GFBC then 
print(`GFBC(G,P,x,K): The conditional generating function (similar to GFB(G,P,x,K)) but conditioned on`):
print(`reaching your goal of having all of the occurrences on a die from the list G taking a total of <=K dice rolls. Try:`):
print(`GFBC([5,3,2],[1/2,1/8,3/8],x,10); `):


elif nargs=1 and args[1]=GFBCgeneral then 
print(`GFBCgeneral(G,P,x,K,g): The conditional generating function (similar to GFBgeneral(G,P,x,K,g)) but conditioned on`):
print(`reaching your goal taking a total of <=K dice rolls where your goal is to have EXACTLY g of the desired occurrences from the list G. Try:`):
print(`GFBCgeneral([5,3,2],[1/2,1/8,3/8],x,10,3);`):


elif nargs=1 and args[1]=GFBCgeneralT then 
print(`GFBCgeneralT(G,P,t,K,g): It computes GFBCgeneral(G,P,x,K,g) and substitutes each variable`):
print(`x[1],x[2],...,x[k] (k is the size of G) to one variable t indicating the total number of dice rolls. Try:`):
print(`GFBCgeneralT([5,3,2],[1/2,1/8,3/8],t,10,3);`):


elif nargs=1 and args[1]=GFBCt then 
print(`GFBCt(G,P,t,K): It computes GFBC(G,P,x,K) and substitutes each variable`):
print(`x[1],x[2],...,x[k] (k is the size of G) to one variable t indicating the total number of dice rolls. Try:`):
print(`GFBCt([5,3,2],[1/2,1/8,3/8],t,10); `):


elif nargs=1 and args[1]=GFBCvg then 
print(`GFBCvg(P,x,S,K): The conditional generating function (similar to GFBvg(P,x,S,K)) but conditioned on`):
print(`reaching your goal taking a total of <=K steps. Try:`):
print(`GFBCvg([1/2,1/8,3/8],x,{{x[1]-5,x[2]-3,x[3]-2}},10); `):


elif nargs=1 and args[1]=GFBvgT then 
print(`GFBvgT(P,t,S,K,x): It computes GFBvg(G,P,x,K,g) and substitutes each variable`):
print(`x[1],x[2],...,x[k] (k is the size of P) to one variable t indicating the total number of occurrences`):
print(`for each x[i]. Try: `):
print(`GFBvgT([1/2,1/2],t,{{2*x[1]+x[2]-7},{x[1]+2*x[2]-7}},5,x); `):


elif nargs=1 and args[1]=GFBgeneral then 
print(`GFBgeneral(G,P,x,K,g): inputs a list G corresponding to the number of (desired) occurrences of each side of a die,`):
print(`a list of probabilities P (for each side of the die), a variable x, an integer K, and an integer g<=nops(G),`):
print(`outputs the prob. generating function in x[1],...,x[k] where k is the number of faces on the die`):
print(`whose coefficient of x[1]^a[1]*...*x[k]^a[k] is the prob. that once you have reached your goal,`):
print(`where your goal is to reach exactly g out of the nops(G) desired ocurrences, the prob.`):
print(`where a[i]=g[i] (when you reached G[i] last) for some i, and and for all j<=g not equal to i, a[j]>=g[j] (when you reach G[j] before reaching G[i]),`):
print(`assuming that you finished in <=K dice rolls and also`):
print(`outputs the probability of not reaching your goal. Try: `):
print(`GFBgeneral([5,3,2],[1/2,1/8,3/8],x,10,3); `):


elif nargs=1 and args[1]=GFBvg then 
print(`GFBvg(P,x,S,K): inputs a probability table P, a variable x, a set S where a member Si of S1`):
print(`is a set of functions in x[1],...,x[k] where k is the size of P, and an integer K`):
print(`outputs the prob. generating function in x[1],...,x[k] whose coefficient of x[1]^a[1]*...*x[k]^a[k]`):
print(`is the prob. that once you have reached your goal,`):
print(`where your goal is to reach a list [b1,b2,...,bk] such that when you plug them into all the sets in S`):
print(`exactly ONE of the sets becomes non-negative, assuming that you finished in <=K steps and also`):
print(`outputs the probability of not reaching your goal. Try: `):
print(`GFBvg([1/2,1/8,3/8],x,{{x[1]-5,x[2]-3,x[3]-1}},10); `):


elif nargs=1 and args[1]=ManyeasySimuG then 
print(`ManyeasySimuG(L,P,K,x): Inputs a list L, a probability table P,`):
print(`a positive integer K and a variable x, outputs a polynomial`):
print(`x[1],...,x[k] where k is the size of L and adds up for each individual`):
print(`the outcome of easySimuG(L,P). Try:`):
print(`ManyeasySimuG([3,4,2],[1/3,1/3,1/3],10,x); `):

 
elif nargs=1 and args[1]=ManygeneralSimuG then 
print(`ManygeneralSimuG(L,P,K,x,g): inputs a list L, a probability table P, a positive integer K,`):
print(`a variable x, and an integer g<=nops(L), outputs a polynomial`):
print(`x[1],...,x[k] where k is the size of L and adds up for each individual`):
print(`the outcome of generalSimuG(L,P,g) and divides by K. Try: `):
print(`ManygeneralSimuG([3,2,2],[1/3,1/3,1/3],10,x,1); `):



elif nargs=1 and args[1]=ManyvgSimuG then 
print(`ManyvgSimuG(P,x,S,K): inputs a probability table P, a variable x,`):
print(`a set S of sets with polynomials, and an integer K, outputs a polynomial`):
print(`x[1],...,x[k] where k is the size of P and adds up for each individual`):
print(`the outcome of vgSimuG(P,x,S,K) and divides by K. Try: `):
print(`ManyvgSimuG([1/3,1/3,1/3],x,{{x[1]+x[2]-10},{x[2]-15},{x[3]-12}},50); `):


elif nargs=1 and args[1]=ManySimuG then 
print(`ManySimuG(L,P,K,x): Inputs a list L, a probability table P,`):
print(`a positive integer K and a variable x, outputs a polynomial`):
print(`x[1],...,x[k] where k is the size of L, and adds up for each individual`):
print(`the outcome of SimuG(L,P). Try: `):
print(`ManySimuG([3,4,2],[1/3,1/3,1/3],10,x); `):


elif nargs=1 and args[1]=OpeAve2 then 
print(`OpeAve2(n,N,p): the linear recurrence operator for the sequence`):
print(`of expected number of coin tosses where the prob. of`):
print(`getting Heads is p and the prob. of getting Tails is 1-p,`):
print(`you finish as soon as the number of Heads is n OR the number of Tails is n. Try:`):
print(`OpeAve2(n,N,p); `):


elif nargs=1 and args[1]=Ope3same then 
print(`Ope3same(n,N): the linear recurrence operator for the sequence`):
print(`of expected number of dice rolls with 3 sides each with prob. 1/3`):
print(`you finish as soon as the number of occurrences for one of the sides equals n. Try:`):
print(`Ope3same(n,N); `):


elif nargs=1 and args[1]=Ope4same then 
print(`Ope4same(n,N): the linear recurrence operator for the sequence`):
print(`of expected number of dice rolls with 4 sides each with prob. 1/4`):
print(`you finish as soon as the number of occurrences for one of the sides equals n. Try:`):
print(`Ope4same(n,N); `):


elif nargs=1 and args[1]=Ope5same then 
print(`Ope5same(n,N): the linear recurrence operator for the sequence`):
print(`of expected number of dice rolls with 5 sides each with prob. 1/5`):
print(`you finish as soon as the number of occurrences for one of the sides equals n. Try:`):
print(`Ope5same(n,N); `):


elif nargs=1 and args[1]=Ope6same then 
print(`Ope6same(n,N): the linear recurrence operator for the sequence`):
print(`of expected number of dice rolls with 6 sides each with prob. 1/6`):
print(`you finish as soon as the number of occurrences for one of the sides equals n. Try:`):
print(`Ope6same(n,N); `):


elif nargs=1 and args[1]=RandomLinear then 
print(`RandomLinear(K1,K2,x,k): Inputs positive integers K1,K2, a symbol x, and a pos. integer k and`):
print(`outputs a linear expression of the form a1x[1]+a2x[2]+....+akx[k]-b`):
print(`where a1, ...,ak are randomly chosen from rand(1..K1) and b from rand(K2..2*K2). Try:`):
print(`RandomLinear(10,50,x,4); `):


elif nargs=1 and args[1]=RandomS then 
print(`RandomS(K1,K2,K3,x,k): Generates random subsets to form a set S to input to GFBCvg(P,x,S,K)`):
print(`where K1 is the range for the coefficients, K2 is the range for b, K3 is the number of subsets`):
print(`and also the number of expressions in each subset, and k is the number of variables. Try: `):
print(`RandomS(10,70,2,x,3); `):



elif nargs=1 and args[1]=Seq2same then 
print(`Seq2same(K,p): Using the second-order recurrence, outputs a list L with the first K values`):
print(`of the expected number of die rolls of a 2-sided die where the prob. of one side is p and the prob. of the other side is 1-p,`):
print(`you finish as soon as the number of occurrences for one of the sides equals i for 1<=i<=K. Try:`):
print(`Seq2same(10,1/2); `):


elif nargs=1 and args[1]=Seq3same then 
print(`Seq3same(K): Using the third-order recurrence, outputs a list L with the first K values`):
print(`of the expected number of die rolls of a 3-sided die where the prob. of each side is 1/3 and`):
print(`you finish as soon as the number of occurrences for one of the sides equals i for 1<=i<=K. Try:`):
print(`Seq3same(10); `):


elif nargs=1 and args[1]=Seq4same then 
print(`Seq4same(K): Using the fourth-order recurrence, outputs a list L with the first K values`):
print(`of the expected number of die rolls of a 4-sided die where the prob. of each side is 1/4 and`):
print(`you finish as soon as the number of occurrences for one of the sides equals i for 1<=i<=K`):
print(`Seq4same(10); `):


elif nargs=1 and args[1]=Seq5same then 
print(`Seq5same(K): Using the fifth-order recurrence, outputs a list L with the first K values`):
print(`of the expected number of die rolls of a 5-sided die where the prob. of each side is 1/5 and`):
print(`you finish as soon as the number of occurrences for one of the sides equals i for 1<=i<=K`):
print(`Seq5same(10); `):


elif nargs=1 and args[1]=Seq6same then 
print(`Seq6same(K): Using the sixth-order recurrence, outputs a list L with the first K values`):
print(`of the expected number of die rolls of a 6-sided die where the prob. of each side is 1/6 and`):
print(`you finish as soon as the number of occurrences for one of the sides equals i for 1<=i<=K`):
print(`Seq6same(10); `):


elif nargs=1 and args[1]=Seqs then 
print(`Seqs(M): the union of all the sequences from Seqs1(M,i) for 1<=i<=nops(M). Try:`):
print(`Seqs([3,1,4]); `):


elif nargs=1 and args[1]=Seqs1 then 
print(`Seqs1(M,i): Given a list M=[m1,m2,...,mk] and a location 1<=i<=k,`):
print(`generates all sequences s of length nops(M) where s[i]=M[i] and s[j]<M[j] j<>i`):
print(`Seqs1([3,1,4],2); `):


elif nargs=1 and args[1]=Simu then 
print(`Simu(p,n,m): inputs a rational number 0<=p<=1, and pos. integers n and m,`):
print(`and with probability p you get Heads and with probility 1-p you get Tails`):
print(`and keeps going until you have reached your goal of having both (at least)`):
print(`n Heads and (at least) m Tails. The output is a pair [i,j] where`):
print(`i is the number of Heads you got and j is the number of Tails you got. Try:`):
print(`Simu(1/2,3,4); `):


elif nargs=1 and args[1]=SimuG then 
print(`SimuG(L,P): Inputs a list of positive integers L, and`):
print(`a list of non-neg. rational numbers P that adds up 1,`):
print(`of the same length as L, and outputs a list X where each values`):
print(`in X is greater or equal to the values in L (pointwise) `):
print(`where X was first a list of zeros and with prob. P[i],`):
print(`X[i] was increased by 1 each time. Try: `):
print(`SimuG([3,4,2], [1/3,1/3,1/3]); `):


elif nargs=1 and args[1]=Stat then 
print(`Stat(F,t,K): Given a generating function F of t, for a r.v. finds the statistical information. Try:`):
print(`Stat((1+t)^10000,t,4); `):


elif nargs=1 and args[1]=StatEasyFS then 
print(`StatEasyFS(G,P,K,L): Inputs a list G of achieving at least G[i] babies of gender i, for AT LEAST ONE OF i=1..k (where k=nops(G)) and the corresponding probability distributions, P,  and a positive integer L`):
print(`outputs`):
print(`[ExpectedFamilySize,StandardDeviation,Skewness, Kurtosis,...] up to the L-th scaled moment for the random variable "total number of babies born" (conditioned on at most K births). Try:`):
print(`StatEasyFS([5,5],[1/2,1/2],6);`):

elif nargs=1 and args[1]=StatFS then 
print(`StatFS(G,P,K,L): Inputs a list G of achieving at least G[i] babies of gender i, for i=1..k (where k=nops(G)) and the corresponding probability distributions, P, a large K (the max. number of births the mother can do), and a positive integer L`):
print(`outputs`):
print(`[ExpectedFamilySize,StandardDeviation,Skewness, Kurtosis,...] up to the L-th scaled moment for the random variable "total number of babies born" (conditioned on at most K births). Try:`):
print(`StatFS([5,5],[1/2,1/2],100,6);`):



elif nargs=1 and args[1]=Trend then 
print(`Trend(p,K,g,K1): inputs a rational number p, 0<p<1, and integers K, g and K1,`):
print(`outputs a list L of length K1 where for each 1<=i<=K1, L[i] is the average family size`):
print(`of having K babies for each gender table G:=[i$d], where d=denom(p), AND the goal is to have `):
print(`EXACTLY g of those genders in G. Try: `):
print(`Trend(1/4,60,4,15); `):


elif nargs=1 and args[1]=vgChop then 
print(`vgChop(d,F,x,S): given a polynomial F in x[1],...,x[d],`):
print(`and a set S where each member Si of S is a set of functions, outputs the part with exponents`):
print(`that satisfy reaching the goal of having one of the sets in S be non-negative. Try:`):
print(`vgChop(3,x[1]^5*x[2]^3*x[3] + x[1]^5*x[2]^3*x[3]^2+x[1]^3*x[2]^2+x[1]^3*x[2]^3,x,{{x[1]-5,x[2]-3,x[3]-1}}); `):


elif nargs=1 and args[1]=vgcomp then 
print(`vgcomp(S,L): Given a set S={S1,S2,...,Sr} such that each Si is a set of functions,`):
print(`and a list L=[L1,L2,...,Lr] both of length r, outputs true if the list L makes all the functions`):
print(`in one of the sets non-negative or false otherwise. Try: `):
print(`vgcomp({{x[1]-3,x[2]-4,x[3]-2}},[3,4,2]); `):


elif nargs=1 and args[1]=vgExpAND then 
print(`vgExpAND(n,P,S,x,K): inputs an integer n, a set S of subsets where each subset,`):
print(`is an expression of the form f(x1,x2,...,xk)=a1x1+a2x2+....+akxk, and a table P with`):
print(`the probabilities such that you get face i with probability Pi`):
print(`and outputs the first n terms of the expected duration of rolling the die.`):
print(`REMARK: If n=1, P=[1/N$N] and S={{x[1],x[2],...,x[N]}} then this reduces to the coupon collector problem. Try:`):
print(`vgExpAND(1,[1/3,1/3,1/3],{{x[1],x[2],x[3]}},x,25);`):


elif nargs=1 and args[1]=vgExpOR then 
print(`vgExpOR(n,P,S,x): inputs an integer n, a set S of subsets where each subset,`):
print(`is an expression of the form f(x1,x2,...,xk)=a1x1+a2x2+....+akxk, and a table P with`):
print(`the probabilities such that you get face i with probability Pi`):
print(`and outputs the FIRST n TERMS of the expected duration of rolling the die. Try:`):
print(`vgExpOR(5,[1/4,1/4,1/4,1/4],{{x[1]},{x[2]},{x[3]},{x[4]}},x); `):


elif nargs=1 and args[1]=vgPGAND then 
print(`vgPGAND(n,P,t,S,x,K): inputs an integer n, a set S of subsets where each subset,`):
print(`is an expression of the form f(x1,x2,...,xk)=a1x1+a2x2+....+akxk, and a table P with`):
print(`the probabilities such that you get face i with probability Pi`):
print(`and outputs the prob. generating function in one variable after substituting`):
print(`each x[i]=t and after subtracting n from each expression (in other words,`):
print(`the goal is to reach the point (y1,y2,....,yk) such that each n<=f(y1,y2,...,yk). Try:`):
print(`vgPGAND(1,[1/2,1/2],t,{{x[1],x[2]}},x,25);`):


elif nargs=1 and args[1]=vgPGOR then 
print(`vgPGOR(n,P,t,S,x): inputs an integer n, a set S of subsets where each subset,`):
print(`is an expression of the form f(x1,x2,...,xk)=a1x1+a2x2+....+akxk, and a table P with`):
print(`the probabilities such that you get face i with probability Pi`):
print(`and outputs the prob. generating function in one variable after substituting`):
print(`each x[i]=t and after subtracting n from each expression (in other words,`):
print(`the goal is to reach the point (y1,y2,....,yk) such that each n<=f(y1,y2,...,yk). Try:`):
print(`vgPGOR(7,[1/2,1/2],t,{{2*x[1]+x[2]},{x[1]+2*x[2]}},x); `):


elif nargs=1 and args[1]=vgSimuG then 
print(`vgSimuG(P,x,S): Inputs a probability table P and an arbitrary set of conditions {S1,S2,...,Sr},`):
print(`where each S[i] is a set of functions of x[1],x[2],...,x[k] such that`):
print(`S[i]={fi1(x[1],x[2],..., x[k]),..., fij(x[1],x[2],..., x[k])} and the size of each S[i]`):
print(`is j=nops(S[i]) and simulates throwing a die with k faces where each face appears with probability P[i]`):
print(`and stops as soon as all the functions in S1 are non-negative OR all the functions in S2 are non-negative, etc`):
print(`(when you plug-in for the symbols x[1],..., x[k] the specific current number of occurrences).`):
print(`When the simulation is done, the output is a list L that contains the number of occurrences`):
print(`where L[i] is the number of times that face i appeared. Try: `):
print(`vgSimuG([1/3,1/3,1/3],x,{{x[1]-3,x[2]-4,x[3]-2}},15); `):

else
print(`There is no Help for`,args):
fi:
 
end:


HelpRec:=proc()
if args=NULL then

print(`HelpRec.txt: A Maple package for empirically guessing partial recurrence`):
print(`equations satisfied by Discrete Functions of ONE Variable`):
print():
print(`For help with a specific procedure, type "HelpRec(procedure_name);"`):
print(`Contains procedures:  `):
print(`findrec, Findrec, FindrecF, SeqFromRec`):
print():

elif nargs=1 and args[1]=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the sequence f of degree DEGREE and order ORDER. Try: `):
print(`findrec([seq(i,i=1..10)],0,2,n,N); `):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC. Try:`):
print(`Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2); `):

elif nargs=1 and args[1]=FindrecF then
print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`):
print(`poly coffs. Try:`):
print(`FindrecF(i->i!,n,N); `):


elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec(N-n-1,n,N,[1],10); `):
fi:
 
end:

#########################
#CheckP(P): Given a probability table outputs true if the entries
# add up to 1
CheckP:=proc(P) local k:
if add(k, k in P)<>1 then
 return false:
fi:
true:
end:


#Simu(p,n,m): inputs a rational number 0<=p<=1, and pos. integers n and m
# and with probability p you get Heads and with probility 1-p you get Tails
# and keeps going until you have reached your goal of having both (at least)
# n Heads and (at least) m Tails. The output is a pair [i,j] where
# i is the number of Heads you got and j is the number of Tails you got
Simu:=proc(p,b,g) local m,n,P,X,k:
if p<0 or p>1 then
 print(`Check that p>=0 and p<=1 `):
 return(FAIL):
fi:
m:=0:
n:=0:
#creates a table with the probabilities
P:=[p,1-p]:
while m<b or n<g do
 X:=RandomVariable(ProbabilityTable(P)):
 k:=convert(Sample(X,1)[1],integer):
 #k will be either 1 or 2 with prob. p and 1-p respectively
 if k=1 then
  m:=m+1:
 else
  n:=n+1:
 fi:
od:
[m,n]:
end:

#comp(L1,L2): inputs two lists L1 and L2 and compares their values
# outputs true if the elements in L2 are greater or equal to
# the elements in L1 (pointwise)
comp:=proc(L1,L2) local n,m,i:
n:=nops(L1):
m:=nops(L2):
if n<>m then
 return(FAIL):
fi:

for i from 1 to n do
 if L2[i]<L1[i] then
  return false:
 fi:
od:
true:
end:


#SimuG(L,P): Inputs a a list of positive integers L, and
# a list of non-neg. rational numbers P that adds up 1, 
# of the same length as L, and outputs a list X where each values
# in X is greater or equal to the values in L (pointwise) 
# where X was first a list of zeros and with prob. P[i],
# X[i] was increased by 1 each time
SimuG:=proc(L,P) local k,k1,p,X,X1:
k1:=nops(P):
if add(p, p in P)<>1 then
 print(`The entries in the probability table should add up to 1`):
 return(FAIL):
fi:

if k1<>nops(L) then
 return(FAIL):
fi:
X:=[0$k1]:

while not comp(L,X) do
 X1:=RandomVariable(ProbabilityTable(P)):
 k:=convert(Sample(X1,1)[1],integer):
 #k will be an integer 1 to k with prob P[i] (1<=i<=k) respectively
 X[k]:=X[k]+1:
od:
X:
end:


#ManySimuG(L,P,K,x): Inputs a list L, a probability table P, 
# a positive integer K and a variable x, outputs a polynomial
# x[1],...,x[k] where k is the size of L, and adds up for each individual
# the outcome of SimuG(L,P)
ManySimuG:=proc(L,P,K,x) local ele,i,k,poly,M:
#L is a list where each element is how many of each gender we want
poly:=0:
for i from 1 to K do
 M:=SimuG(L,P):
 k:=nops(M):
 poly:=poly+mul(x[i]^M[i],i=1..k):
od:
poly/K:
end:


#GFB(G,P,x, K): Inputs a list G corresponding to the number of each desired occurence
# of each side on a die with k:=nops(G) faces, a list of probabilities P (for each side of the die),
# a variable x, and an integer K, outputs the prob. generating function in x[1],...,x[k]
# whose coefficient of x[1]^a[1]*...*x[k]^a[k] is the prob. that once you have reached your goal the prob.
# a[i]=g[i] (when you reached G[i] last) for some i, and for all j not equal to i, a[j]>=g[j] (when you reach G[j] before reaching G[i]),
# assuming that you finished in <=K steps, in other words, you reach all your goals,
# and also outputs the probability of not reaching your goal
GFB:=proc(G,P,x,K) local F,k,AlreadyDone,g1,init,i,j,f1:
k:=nops(G):
if not CheckP(P) or nops(P)<>k then
 print(`Check the probability table or that the size of the probabilities and the list G are the same`):
 return(FAIL):
fi:

if K<add(g1, g1 in G) then
 print(`The sum of the total babies must be greater or equal to the sum of all genders`):
 return(FAIL):
fi: 

init:=add(P[j]*x[j],j=1..k):
F:=1:
AlreadyDone:=0:
for i from 1 to K do
 F:=expand(F*init):
 f1:=Chop(F,x,G):
 AlreadyDone:=expand(AlreadyDone+f1):
 F:=F-f1:
od:
AlreadyDone,subs({seq(x[j]=1,j=1..k)},F):
end:



#Chop(F,x,G): given a polynomial F in x[1],...,x[k] where k is the number of faces on the die,
# and a list G with the number of (desired) occurrences of each side of a die, outputs the part with exponents
# that satisfy reaching the goal of having each of the sides of the die in the list G
Chop:=proc(F,x,G) local k,M,S,m,M1,i:
k:=nops(G):
M:=[op(F)]:
S:=0:
for m in M do
 M1:=[seq(degree(m,x[i]),i=1..k)]:
 if comp(G,M1) then
  S:=S+m:
 fi:
od:
S:
end:


#GFBC(G,P,x,K): The conditional generating function (similar to GFB(G,P,x,K)) but conditioned on
#reaching your goal of having all of the occurrences on a die from the list G taking a total of <=K dice rolls
#this is wanting to reach all your goals
GFBC:=proc(G,P,x,K) local f:
f:=GFB(G,P,x,K):
if f<>FAIL then
 return f[1]/(1-f[2]):
 #1-f[2] is the prob. of reaching your goal
else
 return(FAIL):
fi:
end:


#GFBCt(G,P,t,K): It computes GFBC(G,P,x,K) and substitutes each variable
# x[1],x[2],...,x[k] (k is the size of G) to one variable t indicating the total number of dice rolls
GFBCt:=proc(G,P,t,K) local x,f,k,f1,i:
f:=GFB(G,P,x,K):
if f<>FAIL then
 k:=nops(G):
 f1:=f[1]/(1-f[2]):
 f:=subs({seq(x[i]=t,i=1..k)},f1):
 return f:
fi:
end:


#Stat(F,t,K): Given a generating function F of t, for a r.v. 
# finds the statistical information
Stat:=proc(F,t,K) local gu,mu,sd,k,f:
f:=evalf(F/subs(t=1,F),20):
mu:=subs(t=1,diff(f,t)):
gu:=[mu]:
if K=1 then
 RETURN(gu):
fi:
f:=f/t^mu:
f:=t*diff(t*diff(f,t),t):
sd:=sqrt(subs(t=1,f)):
gu:=[op(gu),sd]:
for k from 3 to K do
 f:=t*diff(f,t):
gu:=[op(gu),subs(t=1,f)/sd^k]:
od:

evalf(gu,20):
end:



########start of easy case

#easycomp(L1,L2): inputs two lists L1 and L2 and compares their values
# outputs true if one of the elements in L2 are greater or equal to
# the corresponding entry in L1 (pointwise)
easycomp:=proc(L1,L2) local n,m,i:
n:=nops(L1):
m:=nops(L2):
if n<>m then
 return(FAIL):
fi:

for i from 1 to n do
 #checks that there is some a[i] in L2 such that b[i]<a[i] for b[i] in L1
 if L2[i]>=L1[i] then
  return true:
 fi:
od:
false:
end:


#easySimuG(L,P): inputs a list of positive integers L, and
# a list of non-neg. rational numbers P that adds up 1, of the same length as L, and
# outputs a list X where there exists one entry in X
# that is greater or equal to the corresponding value in L (pointwise)
# where X was first a list of zeros and with prob. P[i],
# X[i] was increased by 1 each time
easySimuG:=proc(L,P) local k,k1,p,X,X1:
k1:=nops(P):
if add(p, p in P)<>1 then
 print(`The entries in the probability table should add up to 1`):
 return(FAIL):
fi:

if k1<>nops(L) then
 return(FAIL):
fi:
X:=[0$k1]:

while not easycomp(L,X) do
 X1:=RandomVariable(ProbabilityTable(P)):
 k:=convert(Sample(X1,1)[1],integer):
 #k will be an integer 1 to k with prob P[i] (1<=i<=k) respectively
 X[k]:=X[k]+1:
od:
X:
end:


#ManyeasySimuG(L,P,K,x): inputs a list L, a probability table P, 
# a positive integer K and a variable x, outputs a polynomial
# x[1],...,x[k] where k is the size of L and adds up for each individual
# the outcome of easySimuG(L,P)
ManyeasySimuG:=proc(L,P,K,x) local i,k,poly,M:
#L is a list where each element is how many of each gender we want
poly:=0:
for i from 1 to K do
 M:=easySimuG(L,P):
 k:=nops(M):
 poly:=poly+mul(x[i]^M[i],i=1..k):
od:
poly/K:
end:



#easyGFB(G,P,x): inputs a list G corresponding to the number of (desired) occurrences of each side of a die,
# a list of probabilities P (for each side of the die), a variable x,
# outputs the prob. generating function in x[1],...,x[k] where k is the number of faces on the die
# whose coefficient of x[1]^a[1]*...*x[k]^a[k] is the prob. that once you have reached your goal the prob.
# either a[1]=g[1] (if you reached G[1]) or a[2]=g[2] (if you reached G[2]), etc.
#Note: this should be the conditional generating function by default
easyGFB:=proc(G,P,x) local K,g1,F,k,AlreadyDone,g,i,j,f1:
K:=add(g1, g1 in G):
k:=nops(G):
if not CheckP(P) or nops(P)<>k then
 print(`Check the probability table or that the size of the probabilities and the list G are the same`):
 return(FAIL):
fi:

g:=add(P[j]*x[j],j=1..k):
F:=1:
AlreadyDone:=0:
for i from 1 to K do
 F:=expand(F*g):
 f1:=easyChop(F,x,G):
 AlreadyDone:=expand(AlreadyDone+f1):
 F:=F-f1:
od:
AlreadyDone:
end:


#easyChop(F,x,G): given a polynomial F in x[1],...,x[k] where k is the number of faces on the die,
# and a list G with the number of (desired) occurrences of each side of a die, outputs the part with exponents
# that satisfy reaching the goal of reaching all the occurrences of one of the sides of the die from the list G
easyChop:=proc(F,x,G) local k,M,S,m,M1,i:
k:=nops(G):
M:=[op(F)]:
S:=0:
for m in M do
 M1:=[seq(degree(m,x[i]),i=1..k)]:
 if easycomp(G,M1) then
  S:=S+m:
 fi:
od:
S:
end:


#easyGFBt(G,P,t): It computes easyGFB(G,P,x) and substitutes each variable
# x[1],x[2],...,x[k] (k is the size of G) to one variable t indicating the total number of dice rolls
easyGFBt:=proc(G,P,t) local x,f,k,f1,i:
k:=nops(G):
if nops(P)<>k then
  RETURN(FAIL):
fi:
f:=easyGFB(G,P,x):
if f=FAIL then
  RETURN(FAIL):
fi:
expand(subs({seq(x[i]=t,i=1..k)},f)):
end:


#StatEasyFS(G,P,L): Inputs a list G of achieving at least G[i] babies of gender i, for AT LEAST one i (where k=nops(G)) and the corresponding probability distributions, P, and a positive integer L
#outputs
#[ExpectedFamilySize,StandardDeviation,Skewness, Kurtosis,...] up to the L-th scaled moment for the random variable "total number of babies born". Try:
#StatEasyFS([5,5],[1/2,1/2],6);
StatEasyFS:=proc(G,P,L) local t:
Stat(easyGFBt(G,P,t),t,L):
end:


#AllVecs(L): Given a list of positive integers outputs all vectors of positive integers v[i]<=L[i]
AllVecs:=proc(L) local gu,a,L1,gu1:
if nops(L)=1 then
RETURN({seq([a],a=0..L[1])}):
fi:
L1:=[op(1..nops(L)-1,L)]:
gu:=AllVecs(L1):
{seq(seq([op(gu1),a],a=0..L[nops(L)]),gu1 in gu)}:
end:


#easyGFBd(G,P,x): A direct way to compute easyGFB(G,P,x). Try:
#easyGFBd([5,5],[1/2,1/2],x);
easyGFBd:=proc(G,P,x) local i1,i,k,S,v,gu:
k:=nops(G):
if not (nops(G)=nops(P) and convert(P,`+`)=1) then
RETURN(FAIL):
fi:

gu:=0:
for i from 1 to k do
S:=AllVecs([seq(G[i1]-1,i1=1..i-1),seq(G[i1]-1,i1=i+1..k)]):
gu:=gu+
(P[i]*x[i])^G[i]*
add(
mul((P[i1]*x[i1])^v[i1],i1=1..i-1)*mul((P[i1]*x[i1])^v[i1-1],i1=i+1..k)*(convert(v,`+`)+G[i]-1)!/mul(v[i1]!,i1=1..k-1)/(G[i]-1)!,
v in S):
od:
expand(gu):
end:


#order 2
#OpeAve2(n,N,p): the linear recurrence operator for the sequence
# of expected number of coin tosses where the prob. of
# getting Heads is p and the prob. of getting Tails is 1-p,
# you finish as soon as the number of Heads is n OR the number of Tails is n. Try:
# OpeAve2(n,N,p);
OpeAve2:=proc(n,N,p):
-2*(p-1)*(2*n + 1)*p/n + (4*n*p^2 - 4*n*p + 2*p^2 - n - 2*p - 2)*N/(n + 1) + N^2:
end:


#order 3
#Ope3same(n,N): the linear recurrence operator for the sequence
# of expected number of dice rolls with 3 sides each with prob. 1/3
# you finish as soon as the number of occurrences for one of the sides equals n. Try:
# Ope3same(n,N);
Ope3same:=proc(n,N):
-(3*n + 2)*(2*n + 1)*(3*n + 1)*(9*n + 16)/(18*(9*n + 7)*(2 + n)^2*n)
 + (243*n^4 + 1161*n^3 + 2043*n^2 + 1585*n + 458)*N/(9*(n + 1)*(9*n + 7)*(2 + n)^2) 
 - (486*n^3 + 2079*n^2 + 2853*n + 1198)*N^2/(18*(9*n + 7)*(2 + n)^2) + N^3:
end:


#order 4
#Ope4same(n,N): the linear recurrence operator for the sequence
# of expected number of dice rolls with 4 sides each with prob. 1/4
# you finish as soon as the number of occurrences for one of the sides equals n. Try:
# Ope4same(n,N);
Ope4same:=proc(n,N):
(1 + n)*(1 + 2*n)*(3 + 2*n)*(1 + 3*n)*(2 + 3*n)*(1 + 4*n)*(3 + 4*n)*(147600 + 245145*n + 147970*n^2 + 38048*n^3 + 3456*n^4)
-2*n*(3 + 2*n)*(107673150 + 658211935*n + 1778307543*n^2 + 2784317245*n^3 + 2777161071*n^4 + 1821559714*n^5 + 780425472*n^6 + 208857600*n^7 + 31371264*n^8 + 1990656*n^9)*N
+n*(1 + n)*(2185620030 + 11881116237*n + 28149888295*n^2 + 38202277642*n^3 + 32767382228*n^4 + 18431188712*n^5 + 6794048640*n^6 + 1579108608*n^7 + 209129472*n^8 + 11943936*n^9)*N^2
-2*n*(1 + n)*(2 + n)^2*(319483959 + 1357078433*n + 2363003884*n^2 + 2185035148*n^3 + 1162443776*n^4 + 357161472*n^5 + 58761216*n^6 + 3981312*n^7)*N^3
+576*n*(1 + n)*(2 + n)^2*(3 + n)^3*(15833 + 49525*n + 54562*n^2 + 24224*n^3 + 3456*n^4)*N^4:
end:


#order 5
#Ope5same(n,N): the linear recurrence operator for the sequence
# of expected number of dice rolls with 5 sides each with prob. 1/5
# you finish as soon as the number of occurrences for one of the sides equals n. Try:
# Ope5same(n,N);
Ope5same:=proc(n,N):
-((1 + n)*(1 + 2*n)*(3 + 2*n)*(1 + 3*n)*(2 + 3*n)*(1 + 4*n)*(3 + 4*n)*(1 + 5*n)*(2 + 5*n)*(3 + 5*n)*(4 + 5*n)*(-1850089517394854400 -8771409295145570880*n - 18435157669856821680*n^2 - 
22415980500900495912*n^3 - 17103696505089514514*n^4 - 8082097966581118575*n^5 - 1892322845673574750*n^6 + 
249261697425037500*n^7 + 380250989883937500*n^8 + 151378700530546875*n^9 + 34483133651562500*n^10 + 4806119062500000*n^11 + 
383462500000000*n^12 + 13500000000000*n^13)) + n*(3 + 2*n)*(-943343745223441062412800 - 14049607436176624222014720*n - 
96826038622750333708677216*n^2 - 411682252649175960453516240*n^3 - 
1213325329023832975335802752*n^4 - 2637588705548448063937188768*n^5 - 
4388961460077309096668072114*n^6 - 5719703538410196082589470775*n^7 - 
5915207705657944767128748768*n^8 - 4879533424004910472740864672*n^9 - 
3199287206411526483140139600*n^10 - 1640156771294550932197053750*n^11 - 
629640434989234559486217500*n^12 - 158732385850690246647725000*n^13 - 
9823520318963276966718750*n^14 + 12875801061483186533203125*n^15 + 
7318292777116609335937500*n^16 + 2297471437700828906250000*n^17 + 
488974275804132812500000*n^18 + 72428712395781250000000*n^19 + 
7211902781250000000000*n^20 + 436241250000000000000*n^21 + 
12150000000000000000*n^22)*N - 
3*n*(1 + n)*(-8087371549569553121314560 - 112375017735131670031929024*n - 
718385387102793937942799040*n^2 - 2819390701305933800625691584*n^3 - 
7637806149248784102264596856*n^4 - 15205896903280106971714618024*n^5 - 
23099193412785429060054610230*n^6 - 27404472174311014509657531291*n^7 - 
25736250083599258895978854669*n^8 - 19234183385015604052500042097*n^9 - 
11397792014564826109629790225*n^10 - 5263865741053315253630459375*n^11 - 
1808101767986785432621414375*n^12 - 397908771875941272364346875*n^13 - 
11949056145760216320078125*n^14 + 32788259282508692399218750*n^15 + 
16020894212139777812500000*n^16 + 4530951437253154687500000*n^17 + 
880385127100093750000000*n^18 + 119969330569375000000000*n^19 + 
11055330375000000000000*n^20 + 622155000000000000000*n^21 + 
16200000000000000000*n^22)*N^2
+ n*(1 + n)*(2 + n)^2*(-15304920344649197541796608 - 193121282554903079206147296*n - 
1105275248378953218249886008*n^2 - 3827578634564571453568713876*n^3 - 
9011949555082959707257433262*n^4 - 15342399872798386880749372689*n^5 - 
19575518372504432064583588817*n^6 - 19107588694567840522727857019*n^7 - 
14397857314174357358016112325*n^8 - 8352264467286071192731113125*n^9 - 
3654855934228918282968456875*n^10 - 1134765860710751257521978125*n^11 - 
198509303740449027101640625*n^12 + 15033831362055562969531250*n^13 + 
23585328048633443984375000*n^14 + 8687734545716534375000000*n^15 + 
1955779696603718750000000*n^16 + 295460825020625000000000*n^17 + 
29480153625000000000000*n^18 + 1769265000000000000000*n^19 + 
48600000000000000000*n^20)*N^3 
- 2*n*(1 + n)*(2 + n)^2*(3 + n)^3*(-281862096684628550644224 - 3229830082127350390802400*n - 
16476153783696083223375336*n^2 - 49829172921540984717395772*n^3 - 
100133789219361165412231418*n^4 - 141668638331491753561034275*n^5 - 
145451709053975599827811875*n^6 - 109639387557198189496257500*n^7 - 
60239072639752033452475000*n^8 - 23202977555461286665546875*n^9 - 
5472853245092006735546875*n^10 - 250249940446962441406250*n^11 + 
360795049474768359375000*n^12 + 156025596318195312500000*n^13 + 
34724176614531250000000*n^14 + 4654654031250000000000*n^15 + 
357266250000000000000*n^16 + 12150000000000000000*n^17)*N^4
+360000*n*(1 + n)*(2 + n)^2*(3 + n)^3*(4 + n)^4*(-2121206370164352 - 21960005788942200*n - 98699374377269028*n^2 - 
255269554815378606*n^3 - 422755528434011639*n^4 - 469654886584226325*n^5 - 
353753445511399750*n^6 - 175606050109275000*n^7 - 50986508383171875*n^8 - 
3825337547578125*n^9 + 3063348964062500*n^10 + 1257569062500000*n^11 + 
207962500000000*n^12 + 13500000000000*n^13)*N^5:
end:


#order 6
#Ope6same(n,N): the linear recurrence operator for the sequence
# of expected number of dice rolls with 6 sides each with prob. 1/6
# you finish as soon as the number of occurrences for one of the sides equals n. Try:
# Ope6same(n,N);
Ope6same:=proc(n,N):
(1 + n)*(1 + 2*n)*(3 + 2*n)*(5 + 2*n)*(1 + 3*n)*(2 + 3*n)*(4 + 3*n)*(5 + 3*n)*(1 + 4*n)*(3 + 4*n)*(1 + 5*n)*(2 + 5*n)*(3 + 5*n)*(4 + 5*n)*(1 + 6*n)*(5 + 6*n)*(-521906481929159021485167226020864000 - 
 4725009892795393306501288284003655680*n - 
 20589389414212272986480870999804881536*n^2 - 
 57431915510229767764835141467214176992*n^3 - 
 115023108424904492847760165251617103360*n^4 - 
 175831225558051410033420511468090769252*n^5 - 
 212910089043857601013272872892633936436*n^6 - 
 209204662436917314992909796056140214011*n^7 - 
 169551243257919957020869684437452720420*n^8 - 
 114607285738519462424606446795951865886*n^9 - 
 65091426158883583428955613180532444450*n^10 - 
 31204366565466779547405263290929202612*n^11 - 
 12654120942389337087305625187008810770*n^12 - 
 4340855436482707556520966963471018182*n^13 - 
 1256908984742529129115099063331669578*n^14 - 
 305802121010345493685711190037747585*n^15 - 
 62053083868554761985968090923217250*n^16 - 
 10385178076275581256008102239185000*n^17 - 
 1409851003709559171994436572425000*n^18 - 
 151367633870367001091268060000000*n^19 - 12330534874516008177719812500000*
 n^20 - 704604850061840639848125000000*n^21 - 
 22987159364162076187500000000*n^22 - 10864815971512500000000000*n^23 + 
 30587246187300000000000000*n^24 + 860934420000000000000000*n^25)
-n*(3 + 2*n)*(5 + 2*n)*(4 + 3*n)*(5 + 3*n)*(-6508603435344114540075609411253249746739200 - 
 132986121574892939262882018394623955345075200*n - 
 1304853379875731390583118024186707681405292032*n^2 - 
 8223087470524015788154622968115502537147177792*n^3 - 
 37520140897843114287688588393946627108625625536*n^4 - 
 132355738643266288134120672524813958675771873168*n^5 - 
 376066243822530666743919042997195505908956582720*n^6 - 
 884638902297771671388060435245416928243099342764*n^7 - 
 1756277948339832492851720209101048625273972372692*n^8 - 
 2983521325899514741384557454797295761807984855947*n^9 - 
 4380610125036711449513524751551034803840673000106*n^10 - 
 5600478317143071236782481959647850565188038412290*n^11 - 
 6268871196948852074686117263507186814167317418960*n^12 - 
 6168845785625510098349824862731367531873872124377*n^13 - 
 5352750688223570956985586124050012514770037423950*n^14 - 
 4104327254554605632919387199533368019825433580660*n^15 - 
 2784937626453154824233156675535568287233513217660*n^16 - 
 1673518217806245868474656798622983643208235032733*n^17 - 
 890739197511995000622450748874992826047299997542*n^18 - 
 419730132992467134912989272013969543542498684594*n^19 - 
 174904709914276508335285769140212449511082831752*n^20 - 
 64335350186770491186027321566375328921185882975*n^21 - 
 20834000230659728570555722551431407959024002850*n^22 - 
 5918655519103265599215324400109643163567567500*n^23 - 
 1468126189458861217720201449711757002334095000*n^24 - 
 316040250800921323655400598765611715106790000*n^25 - 
 58574625324146060245767175030437320839800000*n^26 - 
 9249643345501638565347578833222287480000000*n^27 - 
 1227102873838057634203873844324644050000000*n^28 - 
 134105310800122595713508265302709000000000*n^29 - 
 11726185641646534961782978955925000000000*n^30 - 
 781980752323613048898827947500000000000*n^31 - 
 36157080564107857717190100000000000000*n^32 - 
 863931543765786640170000000000000000*n^33 + 
 11680069272007770000000000000000000*n^34 + 
 1498858586571024000000000000000000*n^35 + 33473130249600000000000000000000*n^36)*N
+3*n*(1 + n)*(5 + 2*n)*(-3468741505834851266479889502592766732580864000 - 
 67805744473671163269077586399969057722563553280*n - 
 634657661572288025853886393285622371876592756736*n^2 - 
 3808509907846715705707782629263269382650370635648*n^3 - 
 16542890102506193197442847941269219992460435950144*n^4 - 
 55643612503780942124744209050430653773250893149248*n^5 - 
 151309625259650840914367008212044064116069301641680*n^6 - 
 342610921331349381891786634098893742885911951612120*n^7 - 
 659709922292753039089881447097168183429696819206684*n^8 - 
 1096727803404518330401480346090916903135827795766592*n^9 - 
 1591325142743763676801884436528440459914188319458851*n^10 - 
 2030877584262133619560348199133473896285904218656672*n^11 - 
 2292035672009647420453165620465372795604992987905194*n^12 - 
 2296145985490024012755751196707401336363243940925096*n^13 - 
 2047084126314969225553801095994101499184640417722517*n^14 - 
 1626990488446663492676408116323855370916469392464744*n^15 - 
 1154079292593265516354495443897118071880918590531984*n^16 - 
 731079142259128611933700147023056114797085526227040*n^17 - 
 413678635365031011657524650803935322242209601683213*n^18 - 
 209047447995329564742787095357476003940609090919376*n^19 - 
 94283638590477093898437739328942716648331966431794*n^20 - 
 37910478470634768431697257756221225229317849408744*n^21 - 
 13567601920023431499241318606750554713785022787003*n^22 - 
 4312086600279553481892463008528804993565993111440*n^23 - 
 1213423458766971535437712063595788431993296650200*n^24 - 
 301155488041766426441454082680915346208787930000*n^25 - 
 65592691923266110179890720128756991399395050000*n^26 - 
 12457359063688241220652489106065075657446600000*n^27 - 
 2046024603249245617108732798378367330172000000*n^28 - 
 287473744664191323278759646944603786700000000*n^29 - 
 34051822465331049053554260352943824500000000*n^30 - 
 3331333872417339836879740841480760000000000*n^31 - 
 260990140857154881223475450254500000000000*n^32 - 
 15544045644989202403376452650000000000000*n^33 - 
 631161145936019676401851500000000000000*n^34 - 
 11829180924243130481550000000000000000*n^35 + 
 319050902015418510000000000000000000*n^36 + 
 25620984759465360000000000000000000*n^37 + 
 502096953744000000000000000000000*n^38)*N^2
-n*(1 + n)*(2 + n)^2*(-68623702988705174250603876260560400317495296000 - 
 1262669423110430724950567386730637457864380395520*n - 
 11062026217287402034107170544609959014055953054464*n^2 - 
 61815048374136659580968856957093685099471652643840*n^3 - 
 248901864734726760296903546180859476113560228538320*n^4 - 
 773168833197004717241874409191734479342738394335488*n^5 - 
 1936122018219321175934285650243859423869097632613736*n^6 - 
 4029741199661605556135278671616444283474318151702736*n^7 - 
 7126198334922477392310720640253601214681726819596877*n^8 - 
 10878751360936024119050520415938345495069438529706004*n^9 - 
 14499255580830597026313139967413704815489396892287928*n^10 - 
 17003246715078566182355122378428318475073258256415362*n^11 - 
 17634598629138319263183795540233645167812184926373149*n^12 - 
 16227588489565524977650015676204967688440576941169370*n^13 - 
 13275250539453791731584451060123896676992175200474160*n^14 - 
 9665058600243697002927780374857871719146557061798708*n^15 - 
 6265568787331472921686376592178755896265312717183591*n^16 - 
 3616987507531064509284670701621988660465058644149096*n^17 - 
 1858838002417546514561300968149020632309705986673632*n^18 - 
 849888496493071563910200506844359164663551579738154*n^19 - 
 345340484865189665449079412704597928074390234134063*n^20 - 
 124516600894118273180018934042137324971058986548522*n^21 - 
 39753699428906520936682919842919956350953313784080*n^22 - 
 11206420796785796132078301463720928892617892457200*n^23 - 
 2778973131432162806807435700347788373838201030000*n^24 - 
 603312362456455190516425243074992623422825900000*n^25 - 
 113956732771474930255073720348641652377150800000*n^26 - 
 18576458842088062124273061557885026022976000000*n^27 - 
 2585645142312211091914897604059861922100000000*n^28 - 
 302877276449092289080241388347146131000000000*n^29 - 
 29252766822766000982631479809499580000000000*n^30 - 
 2258563060747726952212399202361000000000000*n^31 - 
 132272045888313539897468228700000000000000*n^32 - 
 5258648110716424942432812000000000000000*n^33 - 
 94379772923258600084400000000000000000*n^34 + 
 2765843176279437360000000000000000000*n^35 + 
 209988847613162880000000000000000000*n^36 + 
 4016775629952000000000000000000000*n^37)*N^3
+2*n*(1 + n)*(2 + n)^2*(3 + n)^3*(-2553867699535578137834761688591120783167488000 - 
 43738786419531746973349291608330541459845826560*n - 
 354071651656002313999639017261056283113474287872*n^2 - 
 1816098076735108806521517490284488028294466283136*n^3 - 
 6672304385349166508838143318956794695613353254704*n^4 - 
 18814613800082757823717563783432973889364943462912*n^5 - 
 42588833963645122994594268750478010124425671760024*n^6 - 
 79861435506608663178826858638923858924876237112896*n^7 - 
 126898799159270936366170963176453247390845670503743*n^8 - 
 173644978441350000319785855175003379081838909632246*n^9 - 
 206880644087265088061103017761225559895094869499265*n^10 - 
 216092996916179148794248330998004474339308033358410*n^11 - 
 198667183849552415557823188233478180539884938857682*n^12 - 
 161056652215838855001777937295766720449354931184132*n^13 - 
 115194732239059158940349207862526521917641632237850*n^14 - 
 72674486073088085097684876534827042197820317697428*n^15 - 
 40411854676476472139257154135670983574333429939571*n^16 - 
 19784882860021331086216964341110044553232149113150*n^17 - 
 8515528949859417644930455235729284083415737628989*n^18 - 
 3215855808916955632005749028761438384292041348130*n^19 - 
 1062867447912690095871599545775542248598726578300*n^20 - 
 306413137831721490180752910344569427818398345000*n^21 - 
 76715452732776291476564420010499408486107900000*n^22 - 
 16585195497894191115923707153671413991828300000*n^23 - 
 3072918387249933892315454744713302619716000000*n^24 - 
 483080991520765403249091126408643757100000000*n^25 - 
 63562494773625409777697792821731681000000000*n^26 - 
 6866416217723730695537457964859280000000000*n^27 - 
 591673427376698075676148874882250000000000*n^28 - 
 38770733580938088643650392325000000000000*n^29 - 
 1755625901145017972184679500000000000000*n^30 - 
 40745357038850652313650000000000000000*n^31 + 
 572766199876531890000000000000000000*n^32 + 
 69331499972236080000000000000000000*n^33 + 
 1506290861232000000000000000000000*n^34)*N^4
-144*n*(1 + n)*(2 + n)^2*(3 + n)^3*(4 + n)^4*(-120092897287797746028809268561699170304000 - 
 1916424766057227181041149783548475974041600*n - 
 14334989587302735470575359448172465166466560*n^2 - 
 67418156519003291114003883936229451595784320*n^3 - 
 225534295588091748441352530563064835746359328*n^4 - 
 575572102436994931122995419284118013297936056*n^5 - 
 1173275623591439981091439263163872049034736180*n^6 - 
 1973357637240653091592396449324180404299678294*n^7 - 
 2802866907345483334670906639934579212278581129*n^8 - 
 3415902891637671704474228243408095943259995653*n^9 - 
 3607342212862735395429200102971283824110141720*n^10 - 
 3317499894062110603032354496598977612233790412*n^11 - 
 2660863635620544345593827158187925458331692038*n^12 - 
 1859996739274378062861531601296451974000706166*n^13 - 
 1130929516592202140187177986276783078988408240*n^14 - 
 596568059012382000591655238341247209877041874*n^15 - 
 272177750152264601097703106587988261108128305*n^16 - 
 107018009173233019396613866841656754797135125*n^17 - 
 36107403628547831466919773782938628285797500*n^18 - 
 10397390299989131935922759970511670342572500*n^19 - 
 2537569154946733892089667731071562634325000*n^20 - 
 520097559397661638439037307279488363750000*n^21 - 
 88413827599538717634398515734397012500000*n^22 - 
 12250797350315638869425505929796000000000*n^23 - 
 1348699899654438721450450643981250000000*n^24 - 
 113270290235215191093912249375000000000*n^25 - 
 6735452699757109546100025000000000000*n^26 - 
 235595522897457000142500000000000000*n^27 - 
 1060846498862365500000000000000000*n^28 + 
 274295255893956000000000000000000*n^29 + 8368282562400000000000000000000*n^30)*N^5
+233280000*n*(1 + n)*(2 + n)^2*(3 + n)^3*(4 + n)^4*(5 + n)^5*(-7487404462478851637438947737600 - 
 111123577892325785537668487915520*n - 765394449369850160216790048149760*n^2 - 3283693320917400108869063817301056*n^3 - 
 9934772477890975548560841748833632*n^4 - 
 22759319988067243820131118737213320*n^5 - 
 41394651865241091785181959392512428*n^6 - 
 61824862486964440022412619925278758*n^7 - 
 77640888994400157990397740633268391*n^8 - 
 83206977327277980774103472605873545*n^9 - 
 76628812557395358014028482523784309*n^10 - 
 60687822250116890364750785968030701*n^11 - 
 41189825729281006908544829663669827*n^12 - 
 23819461482101139146053735682380515*n^13 - 
 11653019864085868594485343346225803*n^14 - 
 4783776054270853660773079256631585*n^15 - 
 1632200463073651014638831735597250*n^16 - 
 457383790999270190821399070535000*n^17 - 103625546492467986428164057425000*n^18 - 18567546243814516963728060000000*n^19 - 
 2545346615001595902721687500000*n^20 - 252653908361162826223125000000*n^21 - 16275337815122488687500000000*n^22 - 486678398466712500000000000*n^23 + 9063885687300000000000000*n^24
 + 860934420000000000000000*n^25)*N^6:
end:


#Seq2same(K,p): Using the second-order recurrence, outputs a list L with the first K values 
# of the expected number of die rolls of a 2-sided die where the prob. of one side is p and the prob. of the other side is 1-p,
# you finish as soon as the number of occurrences for one of the sides equals i for 1<=i<=K
Seq2same:=proc(K,p) local n,N,INI,ope:
INI:=[1, -2*p^2 + 2*p + 2]:
ope:=OpeAve2(n,N,p):
SeqFromRec(ope,n,N,INI,K):
end:


#Seq3same(K): Using the third-order recurrence, outputs a list L with the first K values
# of the expected number of die rolls of a 3-sided die where the prob. of each side is 1/3 and
# you finish as soon as the number of occurrences for one of the sides equals i for 1<=i<=K
Seq3same:=proc(K) local n,N,INI,ope:
INI:=[1, 26/9, 409/81]:
ope:=Ope3same(n,N):
SeqFromRec(ope,n,N,INI,K):
end:


#Seq4same(K): Using the fourth-order recurrence, outputs a list L with the first K values
# of the expected number of die rolls of a 4-sided die where the prob. of each side is 1/4 and
# you finish as soon as the number of occurrences for one of the sides equals i for 1<=i<=K
Seq4same:=proc(K) local n,N,INI,ope:
INI:=[1, 103/32, 48039/8192, 2288659/262144]:
ope:=Ope4same(n,N):
SeqFromRec(ope,n,N,INI,K):
end:


#Seq5same(K): Using the fifth-order recurrence, outputs a list L with the first K values
# of the expected number of die rolls of a 5-sided die where the prob. of each side is 1/5 and
# you finish as soon as the number of occurrences for one of the sides equals i for 1<=i<=K
Seq5same:=proc(K) local n,N,INI,ope:
INI:=[1, 2194/625, 2580591/390625, 2443849764/244140625, 2076630126481/152587890625]:
ope:=Ope5same(n,N):
SeqFromRec(ope,n,N,INI,K):
end:


#Seq6same(K): Using the sixth-order recurrence, outputs a list L with the first K values
# of the expected number of die rolls of a 6-sided die where the prob. of each side is 1/6 and
# you finish as soon as the number of occurrences for one of the sides equals i for 1<=i<=K
Seq6same:=proc(K) local n,N,INI,ope:
INI:=[1, 1223/324, 4084571/559872, 247150321423/22039921152, 56252877655712005/3656158440062976, 2597868106693535971/131621703842267136]:
ope:=Ope6same(n,N):
SeqFromRec(ope,n,N,INI,K):
end:


#Ck(k,K): inputs integers k and K where 3<=k<=6 and outputs the constant Ck where
# Ck is asymptotically equal to (k*N-ak(K))/(k*sqrt(K)) as K goes to infinity
# where ak(K) is the expected number of die rolls of a k-sided die such that
# your goal is to see K ocurrences of one side and each side has prob. 1/k
Ck:=proc(k,K) local a:
if k<3 or k>6 then
 print(`the first input can only be 3,4,5 or 6`):
 return(FAIL):
elif k=3 then
 a:=Seq3same(K):
 return (k*K-a[K])/(k*sqrt(K)):
elif k=4 then
 a:=Seq4same(K):
 return (k*K-a[K])/(k*sqrt(K)):
elif k=5 then
 a:=Seq5same(K):
 return (k*K-a[K])/(k*sqrt(K)):
elif k=6 then
 a:=Seq6same(K):
 return (k*K-a[K])/(k*sqrt(K)):
fi:
end:

#######end of easy case


#StatFS(G,P,K,L): Inputs a list G of achieving at least G[i] babies of gender i, for i=1..k (where k=nops(G)) and the corresponding probability distributions, P, a large K (the max. number of births the mother can do), and a positive integer L
#outputs
#[ExpectedFamilySize,StandardDeviation,Skewness, Kurtosis,...] up to the L-th scaled moment for the random variable "total number of babies born" (conditioned on at most K births). Try:
#StatFS([5,5],[1/2,1/2],100,6);
StatFS:=proc(G,P,K,L) local t:
Stat(GFBCt(G,P,K,t),t,L):
end:


#CheckEasyGenderRatioConjecture(G,P): Inputs a list G corresponding to the number of each genders, and
# a list of probabilities P and outputs the ratio of the expected babies from each gender
# and their corresponding probability such that there is only one of the genders that was met
CheckEasyGenderRatioConjecture:=proc(G,P) local K,k,f,x,i,L1,f1,j,g1:
K:=add(g1, g1 in G):
k:=nops(G):
if not CheckP(P) or nops(P)<>k then
 print(`Check the probability table or that the size of the probabilities and the list G are the same`):
 return(FAIL):
fi:
f:=easyGFB(G,P,x):
L1:=Array([0$k]): 

if f<>FAIL and subs({seq(x[i]=1,i=1..k)},f)=1 then
 for i from 1 to k do
  f1:=diff(f,x[i]):
  L1[i]:=subs({seq(x[j]=1,j=1..k)},f1):
 od:
 L1:=convert(L1,list):
 #outputs the ratio of the probability and the expected babies from each gender
 return {seq(P[i]/L1[i],i=1..k)}:
fi:
end:


########start of general 
#generalcomp(L1,L2,g): inputs two lists L1 and L2 (same length) and an integer g<=nops(L1)
# compares their values and outputs true if there are g elements in L1 and L2 such that
# those elements in L2 are greater or equal to the corresponding entry in L1 (pointwise)
generalcomp:=proc(L1,L2,g) local n,m,i,co:
n:=nops(L1):
m:=nops(L2):
if n<>m or g>m then
 return(FAIL):
fi:
co:=0:
for i from 1 to n do
 #checks that there is some a[i] in L2 such that b[i]<a[i] for b[i] in L1
 if L2[i]>=L1[i] then
  co:=co+1:
 fi:
od:

if co=g then
 return true:
fi:
false:
end:


#generalSimuG(L,P,g): inputs a list of positive integers L,
# a list of non-neg. rational numbers P that adds up 1, of the same length as L, and
# an integer g<=nops(L) and outputs a list X where there exists EXACTLY g entries in X
# that are greater or equal to the corresponding value in L (pointwise)
# where X was first a list of zeros and with prob. P[i],
# X[i] was increased by 1 each time
generalSimuG:=proc(L,P,g) local k,k1,p,X,X1:
k1:=nops(P):
if add(p, p in P)<>1 then
 print(`The entries in the probability table should add up to 1`):
 return(FAIL):
fi:

if k1<>nops(L) then
 return(FAIL):
fi:
X:=[0$k1]:

while not generalcomp(L,X,g) do
 X1:=RandomVariable(ProbabilityTable(P)):
 k:=convert(Sample(X1,1)[1],integer):
 #k will be an integer 1 to k with prob P[i] (1<=i<=k) respectively
 X[k]:=X[k]+1:
od:
X:
end:


#ManygeneralSimuG(L,P,K,x,g): inputs a list L, a probability table P, a positive integer K,
# a variable x, and an integer g<=nops(L), outputs a polynomial
# x[1],...,x[k] where k is the size of L and adds up for each individual
# the outcome of generalSimuG(L,P,g) and divides by K
ManygeneralSimuG:=proc(L,P,K,x,g) local ele,i,k,poly,M:
#L is a list where each element is how many of each gender we want
poly:=0:
for i from 1 to K do
 M:=generalSimuG(L,P,g):
 k:=nops(M):
 poly:=poly+mul(x[i]^M[i],i=1..k):
od:
poly/K:
end:


#GFBgeneral(G,P,x,K,g): inputs a list G corresponding to the number of (desired) occurrences of each side of a die,
# a list of probabilities P (for each side of the die), a variable x, an integer K, and an integer g<=nops(G),
# outputs the prob. generating function in x[1],...,x[k] where k is the number of faces on the die
# whose coefficient of x[1]^a[1]*...*x[k]^a[k] is the prob. that once you have reached your goal,
# where your goal is to reach exactly g out of the nops(G) desired ocurrences, the prob.
# where a[i]=g[i] (when you reached G[i] last) for some i, and and for all j<=g not equal to i, a[j]>=g[j] (when you reach G[j] before reaching G[i]),
# assuming that you finished in <=K dice rolls and also
# outputs the probability of not reaching your goal
GFBgeneral:=proc(G,P,x,K,g) local F,k,AlreadyDone,g1,init,i,j,f1:
k:=nops(G):
if not CheckP(P) or nops(P)<>k or k<g then
 print(`Check the probability table or that the size of the probabilities and the list G are the same`):
 print(`or that the last input is less than or equal to the size of G`):
 return(FAIL):
fi:

if K<add(g1, g1 in G) then
 print(`The sum of the total babies must be greater or equal to the sum of all genders`):
 return(FAIL):
fi: 

init:=add(P[j]*x[j],j=1..k):
F:=1:
AlreadyDone:=0:
for i from 1 to K do
 F:=expand(F*init): #this is P(x1,x2,...,xk)S_{K-1}(x1,x2,...,xk)
 f1:=generalChop(F,x,G,g): #this is N_K(x1,x2,...,xk) in the paper
 AlreadyDone:=expand(AlreadyDone+f1): #this is F_K(x1,x2,...,xk)=F_{K-1}+N_K(x1,x2,...,xk)
 F:=F-f1: #this is S_K(x1,x2,...,xk)
od:
AlreadyDone, subs({seq(x[j]=1,j=1..k)},F):
end:



#generalChop(F,x,G,g): given a polynomial F in x[1],...,x[k] where k is the number of faces on the die,
# and a list G with number of (desired) occurrences of each side of a die, outputs the part with exponents
# that satisfy reaching the goal of having exactly g out of the nops(G) desired ocurrences
generalChop:=proc(F,x,G,g) local k,M,S,m,M1,i:
k:=nops(G):
M:=[op(F)]:
S:=0:
for m in M do
 M1:=[seq(degree(m,x[i]),i=1..k)]:
 if generalcomp(G,M1,g) then
  S:=S+m:
 fi:
od:
S:
end:


#GFBCgeneral(G,P,x,K,g): The conditional generating function (similar to GFBgeneral(G,P,x,K,g)) but conditioned on
#reaching your goal taking a total of <=K dice rolls
GFBCgeneral:=proc(G,P,x,K,g) local f:
f:=GFBgeneral(G,P,x,K,g):
if f<>FAIL then
 return f[1]/(1-f[2]):
 #1-f[2] is the prob. of reaching your goal
else
 return(FAIL):
fi:
end:


#GFBCgeneralT(G,P,t,K,g): It computes GFBCgeneral(G,P,x,K,g) and substitutes each variable
# x[1],x[2],...,x[k] (k is the size of G) to one variable t indicating the total number of dice rolls
GFBCgeneralT:=proc(G,P,t,K,g) local x,f,k,f1,i:
f:=GFBgeneral(G,P,x,K,g):
if f<>FAIL then
 k:=nops(G):
 f1:=f[1]/(1-f[2]):
 f:=subs({seq(x[i]=t,i=1..k)},f1):
 return f:
fi:
end:


#Trend(p,K,g,K1): inputs a rational number p, 0<p<1, and integers K, g and K1,
# outputs a list L of length K1 where for each 1<=i<=K1, L[i] is the average family size
# of having K babies for each gender table G:=[i$d], where d=denom(p), AND the goal is to have 
# EXACTLY g of those genders in G 
Trend:=proc(p,K,g,K1) local L,L1,L2,P,i,k,x:
k:=denom(p):
if K1<k*K1 then
 print(`Check that K (second input) is greater or equal to`, k*K1):
 return FAIL:
fi:
P:=[p$k]:
L:=[seq(GFBCgeneralT([i$k],P,x,K,g),i=1..K1)]:
L1:=[seq(subs(x=1,diff(L[i],x)),i=1..K1)]:
L2:=evalf([seq(L1[i]/i,i=1..K1)]);
end:


#CheckGenderRatioConjecture(G,P,K): Inputs a list G corresponding to the number of each genders,
# a list of probabilities P and an integer K and outputs the ratio of the expected babies from each gender
# and their corresponding probability for the general case in which every goal is met
CheckGenderRatioConjecture:=proc(G, P, K) local k,f,x,i,L1,f1,j:
k:=nops(G):
if not CheckP(P) or nops(P)<>k then
 print(`Check the probability table or that the size of the probabilities and the list G are the same`):
 return(FAIL):
fi:
f:=GFBC(G,P,x,K):
L1:=Array([0$k]):
Digits:=20:
if f<>FAIL and subs({seq(x[i]=1,i=1..k)},f)=1 then
 for i from 1 to k do
  f1:=diff(f,x[i]):
  L1[i]:=evalf(subs({seq(x[j]=1,j=1..k)},f1)):
 od:
 L1:=convert(L1,list):
 #outputs the ratio of the probability and the expected babies from each gender
 return {seq(evalf(P[i]/L1[i]),i=1..k)}:
fi:
end:
########end of general


########start of super general
#CheckGenderRatioConjectureVG(S,P,K): Inputs a set S where a member Si of S1
# is a set of functions in x[1],...,x[k] and k is the size of P
# a list of probabilities P and an integer K and outputs the ratio of the expected babies from each gender
# and their corresponding probability for the general case in which exactly one set in S is non-negative
CheckGenderRatioConjectureVG:=proc(S, P, K) local k,f,x,i,L1,f1,j:
k:=nops(P):
if not CheckP(P) or nops(P)<>k then
 print(`Check the probability table or that the size of the probabilities and the list G are the same`):
 return(FAIL):
fi:
f:=GFBCvg(P,x,S,K):
L1:=Array([0$k]):
Digits:=20:
if f<>FAIL and subs({seq(x[i]=1,i=1..k)},f)=1 then
 for i from 1 to k do
  f1:=diff(f,x[i]):
  L1[i]:=evalf(subs({seq(x[j]=1,j=1..k)},f1)):
 od:
 L1:=convert(L1,list):
 #outputs the ratio of the probability and the expected babies from each gender
 return {seq(evalf(P[i]/L1[i]),i=1..k)}:
fi:
end:


#RandomLinear(K1,K2,x,k): Inputs positive integers K1,K2, a symbol x, and a pos. integer k and
# outputs a linear expression of the form a1x[1]+a2x[2]+....+akx[k]-b
# where a1, ...,ak are randomly chosen from rand(1..K1) and b from rand(K2..2*K2)
RandomLinear:=proc(K1,K2,x,k) local L,r1,b,i:
L:=[]:
r1:=rand(1..K1):
b:=rand(K2..2*K2)():
for i from 1 to k do
 L:=[op(L), r1()]:
od:
add(L[i]*x[i],i=1..k)-b:
end:


#RandomS(K1,K2,K3,x,k): Generates random subsets to form a set S to input to GFBCvg(P,x,S,K)
# where K1 is the range for the coefficients, K2 is the range for b, K3 is the number of subsets
# and also the number of expressions in each subset, and k is the number of variables.
RandomS:=proc(K1,K2,K3,x,k) local poly,S,k1,k2,i,i1:
S:={}:
for i from 1 to K3 do
 S:=S union {{seq(RandomLinear(K1,K2,x,k),i1=1..K3)}}:
od:
S:
end:


#vgcomp(S,L): Given a set S={S1,S2,...,Sr} such that each Si is a set of functions,
# and a list L=[L1,L2,...,Lr] both of length r, outputs true if the list L makes all the functions
# in one of the sets non-negative or false otherwise
vgcomp:=proc(S,L) local i,j,j1,i1,k,k1,L1,ele:
k:=nops(L):
k1:=nops(S):
L1:=[]:
for i from 1 to k1 do
 j:=nops(S[i]):
 L1:=[op(L1),[seq(subs({seq(x[i1]=L[i1],i1=1..k)},S[i][j1]),j1=1..j)]]:
od:

for ele in L1 do
  if generalcomp([seq(0,i1=1..nops(ele))],ele,nops(ele)) then
   return true:
  fi:
od:
false:
end:


#vgSimuG(P,x,S): Inputs an arbitrary set of conditions {S1,S2,...,Sr},
# where each S[i] is a set of functions of x[1],x[2],...,x[k] such that
# S[i]={fi1(x[1],x[2],..., x[k]),..., fij(x[1],x[2],..., x[k])} and the size of each S[i]
# is j=nops(S[i]) and simulates throwing a die with k faces and stops as soon as
# (when you plug-in for the symbols x[1],..., x[k] the specific current number of occurrences)
# all the functions in S1 are non-negative OR all the functions in S2 are non-negative, etc
# returns a list L of the number of occurrences where L[i] is the number of times that face i appeared
vgSimuG:=proc(P,x,S,K) local k,k1,X,X1,s:
k1:=nops(P):
X:=[0$k1]:
s:=0:
while not vgcomp(S,X) and s<=K do
 X1:=RandomVariable(ProbabilityTable(P)):
 k:=convert(Sample(X1,1)[1],integer):
 #k will be an integer 1 to k with prob P[i] (1<=i<=k) respectively
 X[k]:=X[k]+1:
 s:=s+1:
od:

if s<=K and vgcomp(S,X) then
 return(X):
else
 return(FAIL):
fi:
end:


#ManyvgSimuG(P,x,S,K): inputs a probability table P, a variable x,
# a set S of sets with polynomials, and an integer K, outputs a polynomial
# x[1],...,x[k] where k is the size of P and adds up for each individual
# the outcome of vgSimuG(P,x,S,K) and divides by K
ManyvgSimuG:=proc(P,x,S,K) local ele,i,k,poly,M:
#L is a list where each element is how many of each gender we want
poly:=0:
for i from 1 to K do
 M:=vgSimuG(P,x,S,K):
 if M<>FAIL then
  k:=nops(M):
  poly:=poly+mul(x[i]^M[i],i=1..k):
 fi:
od:
poly/K:
end:


#GFBvg(P,x,S,K): inputs a probability table P, a variable x, a set S where a member Si of S1
# is a set of functions in x[1],...,x[k] where k is the size of P, and an integer K-1
# outputs the prob. generating function in x[1],...,x[k] whose coefficient of x[1]^a[1]*...*x[k]^a[k]
# is the prob. that once you have reached your goal, 
# where your goal is to reach a list [b1,b2,...,bk] such that when you plug them into all the sets in S1
# exactly ONE of the sets becomes non-negative, assuming that you finished in <=K steps and also
# outputs the probability of not reaching your goal
GFBvg:=proc(P,x,S,K) local k1,init,j,F,AlreadyDone,f1,i:
k1:=nops(P):
if not CheckP(P) or nops(P)<>k1 then
 print(`Check the probability table or that the size of the probabilities and that the number of variables are the same`):
 return(FAIL):
fi:

init:=add(P[j]*x[j],j=1..k1):
F:=1:
AlreadyDone:=0:
for i from 1 to K do
 F:=expand(F*init): #this is S_{K-1}(x1,x2,...,xk)*P(x1,x2,...,xk)
 f1:=vgChop(k1,F,x,S): #this is N_K(x1,x2,...,xk)
 AlreadyDone:=expand(AlreadyDone+f1): #this is G_K(x1,x2,...,xk)=G_{K-1}+N_K(x1,x2,...,xk)
 F:=F-f1: #this is S_K(x1,x2,...,xk)
od:
AlreadyDone, subs({seq(x[j]=1,j=1..k1)},F):
end:


#GFBCvg(P,x,S,K): The conditional generating function (similar to GFBvg(P,x,S,K)) but conditioned on
# reaching your goal taking a total of <=K steps
GFBCvg:=proc(P,x,S,K) local f:
f:=GFBvg(P,x,S,K):
if f[2]<>1 then
 return f[1]/(1-f[2]):
 #1-f[2] is the prob. of reaching your goal
else
 print(`Increase K (the last entry)`):
 return(FAIL):
fi:
end:


#d will be the number of sides in the die
#vgChop(d,F,x,S): given a polynomial F in x[1],...,x[d],
# and a set S where each member Si of S is a set of functions, outputs the part with exponents
# that satisfy reaching the goal of having one of the sets in S be non-negative
vgChop:=proc(d,F,x,S) local M,S1,m,M1,i:
M:=[op(F)]:
S1:=0:
for m in M do
 M1:=[seq(degree(m,x[i]),i=1..d)]:
 if vgcomp(S,M1) then
  S1:=S1+m:
 fi:
od:
S1:
end:


#GFBvgT(P,t,S,K,x): It computes GFBvg(G,P,x,K,g) and substitutes each variable
# x[1],x[2],...,x[k] (k is the size of P) to one variable t indicating the total number of occurrences
# for each x[i]. Try:
# GFBvgT([1/2,1/2],t,{{2*x[1]+x[2]-7},{x[1]+2*x[2]-7}},5,x);
GFBvgT:=proc(P,t,S,K,x) local f,k,f1,i:
f:=GFBvg(P,x,S,K)[1]:
if f<>FAIL then
 k:=nops(P):
 f1:=subs({seq(x[i]=t,i=1..k)},f):
 return f1:
fi:
end:


#vgPGOR(n,P,t,S,x): inputs an integer n, a set S of subsets where each subset,
# is an expression of the form f(x1,x2,...,xk)=a1x1+a2x2+....+akxk, and a table P with
# the probabilities such that you get face i with probability Pi
# and outputs the prob. generating function in one variable after substituting
# each x[i]=t and after subtracting n from each expression (in other words,
# the goal is to reach the point (y1,y2,....,yk) such that each n<=f(y1,y2,...,yk). Try:
# vgPGOR(7,[1/2,1/2],t,{{2*x[1]+x[2]},{x[1]+2*x[2]}},x);
vgPGOR:=proc(n,P,t,S,x) local S2,s,eq,K,i:
S2:={seq( { seq(eq - n, eq in s) }, s in S )};
for i from 1 while GFBvgT(P,1,S2,i,x)<>1 do
od:
K:=i:
GFBvgT(P,t,S2,K,x):
end:


#vgExpOR(n,P,S,x): inputs an integer n, a set S of subsets where each subset,
# is an expression of the form f(x1,x2,...,xk)=a1x1+a2x2+....+akxk, and a table P with
# the probabilities such that you get face i with probability Pi
# and outputs the first n terms of the expected duration of rolling the die. Try:
# vgExp(20,[1/2,1/2],{{2*x[1]+x[2]},{x[1]+2*x[2]}});
vgExpOR:=proc(n,P,S,x) local L,len,i,j,i1,K,t:
L:=[seq(vgPGOR(j,P,t,S,x),j=1..n)]:
len:=nops(L):
[seq(subs(t=1,diff(L[i1],t)),i1=1..len)]
end:


#vgPGAND(n,P,t,S,x,K): inputs an integer n, a set S of subsets where each subset,
# is an expression of the form f(x1,x2,...,xk)=a1x1+a2x2+....+akxk, and a table P with
# the probabilities such that you get face i with probability Pi
# and outputs the prob. generating function in one variable after substituting
# each x[i]=t and after subtracting n from each expression (in other words,
# the goal is to reach the point (y1,y2,....,yk) such that each n<=f(y1,y2,...,yk). Try:
# vgPGAND(1,[1/2,1/2],t,{{x[1],x[2]}},x,100);
vgPGAND:=proc(n,P,t,S,x,K) local S2,s,eq,i:
S2:={seq( { seq(eq - n, eq in s) }, s in S )};
GFBvgT(P,t,S2,K,x):
end:


#vgExpAND(n,P,S,x,K): inputs an integer n, a set S of subsets where each subset,
# is an expression of the form f(x1,x2,...,xk)=a1x1+a2x2+....+akxk, and a table P with
# the probabilities such that you get face i with probability Pi
# and outputs the first n terms of the expected duration of rolling the die.
# REMARK: If n=1, P=[1/N$N] and S={{x[1],x[2],...,x[N]}} then this reduces to the coupon collector problem. Try:
# vgExpAND(1,[1/3,1/3,1/3],{{x[1],x[2],x[3]}},x,25); 
vgExpAND:=proc(n,P,S,x,K) local L,len,j,i1,t:
L:=[seq(vgPGAND(j,P,t,S,x,K),j=1..n)]:
len:=nops(L):
[seq(subs(t=1,diff(L[i1],t)),i1=1..len)]
end:


########end of super general


######## (faster) special cases:


F1:=proc(m1,m2,m3,p1,p2,p3,t) local x1,x2,x3:
add(add((x1+x2+m3-1)!/x1!/x2!/(m3-1)!*p1^x1*p2^x2*p3^m3*t^(x1+x2+m3),x1=0..m1-1),x2=0..m2-1):
end:

F2:=proc(m1,m2,m3,p1,p2,p3,t) local x1,x2,x3:
add(add((x1+m2-1+x3)!/x1!/x3!/(m2-1)!*p1^x1*p2^m2*p3^x3*t^(x1+m2+x3),x1=0..m1-1),x3=0..m3-1):
end:

F3:=proc(m1,m2,m3,p1,p2,p3,t) local x1,x2,x3:
add(add((m1-1+x2+x3)!/(m1-1)!/x2!/x3!*p1^m1*p2^x2*p3^x3*t^(m1+x2+x3),x2=0..m2-1),x3=0..m3-1):
end:


#F(m1,m2,m3,p1,p2,p3,t): Inputs non-negative integers m1, m2, and m3,
# probabilities p1,p2,p3, and a variable t, outputs the 
# probability generating function of the duration of the game in which
# with prob. p1 the die lands on side 1, with prob. p2 the die lands on side 2,
# with prob. p3 the die lands on side 3 and 
# at each die roll, you keep track of the number of appearances of each side,
# the game stops until either the number of appearances of side 1 equals m1,
# OR the number of appearances of side 2 equals m2, OR
# the number of appearances of side 3 equals m3  
F:=proc(m1,m2,m3,p1,p2,p3,t):
F1(m1,m2,m3,p1,p2,p3,t)+F2(m1,m2,m3,p1,p2,p3,t)+F3(m1,m2,m3,p1,p2,p3,t):
end:



#ExpF(n,p1,p2,p3): outputs the first k terms of the expected duration of the game
# in which with prob. p1 the die lands on side 1, with prob. p2 the die lands on side 2,
# with prob. p3 the die lands on side 3 and
# at each die roll, you keep track of the number of appearances of each side,
# the game stops until either the number of appearances of side 1 equals k,
# OR the number of appearances of side 2 equals k, OR
# the number of appearances of side 3 equals k, where 1<=k<=n  
ExpF:=proc(n,p1,p2,p3) local t,k:
[seq(subs(t=1,diff(F(k,k,k,p1,p2,p3,t),t)),k=1..n)]
end:

G1:=proc(m1,m2,m3,m4,p1,p2,p3,p4,t) local x1,x2,x3,x4:
add(add(add((x1+x2+x3+m4-1)!/x1!/x2!/x3!/(m4-1)!*p1^x1*p2^x2*p3^x3*p4^m4*t^(x1+x2+x3+m4),x1=0..m1-1),x2=0..m2-1),x3=0..m3-1):
end:

G2:=proc(m1,m2,m3,m4,p1,p2,p3,p4,t) local x1,x2,x3,x4:
add(add(add((x1+x2+m3-1+x4)!/x1!/x2!/x4!/(m3-1)!*p1^x1*p2^x2*p3^m3*p4^x4*t^(x1+x2+m3+x4),x1=0..m1-1),x2=0..m2-1),x4=0..m4-1):
end:

G3:=proc(m1,m2,m3,m4,p1,p2,p3,p4,t) local x1,x2,x3,x4:
add(add(add((x1+m2-1+x3+x4)!/x1!/x3!/x4!/(m2-1)!*p1^x1*p2^m2*p3^x3*p4^x4*t^(x1+m2+x3+x4),x1=0..m1-1),x3=0..m3-1),x4=0..m4-1):
end:

G4:=proc(m1,m2,m3,m4,p1,p2,p3,p4,t) local x1,x2,x3,x4:
add(add(add((m1-1+x2+x3+x4)!/x2!/x3!/x4!/(m1-1)!*p1^m1*p2^x2*p3^x3*p4^x4*t^(m1+x2+x3+x4),x2=0..m2-1),x3=0..m3-1),x4=0..m4-1):
end:


#G(m1,m2,m3,m4,p1,p2,p2,p4,t): Inputs non-negative integers m1, m2, m3, and m4
# probabilities p1,p2,p3,p4, and a variable t, outputs the 
# probability generating function of the duration of the game in which
# with prob. p1 the die lands on side 1, with prob. p2 the die lands on side 2,
# with prob. p3 the die lands on side 3, with prob. p4 the die lands on side 4,
# and at each die roll, you keep track of the number of appearances of each side,
# the game stops until either the number of appearances of side 1 equals m1,
# OR the number of appearances of side 2 equals m2, OR
# the number of appearances of side 3 equals m3 OR
# the number of appearances of side 4 equals m4
G:=proc(m1,m2,m3,m4,p1,p2,p3,p4,t):
G1(m1,m2,m3,m4,p1,p2,p3,p4,t)+G2(m1,m2,m3,m4,p1,p2,p3,p4,t)+G3(m1,m2,m3,m4,p1,p2,p3,p4,t)+G4(m1,m2,m3,m4,p1,p2,p3,p4,t):
end:



#ExpG(n,p1,p2,p3,p4): outputs the first n terms of the expected duration of the game
# in which with prob. p1 the die lands on side 1, with prob. p2 the die lands on side 2,
# with prob. p3 the die lands on side 3, with prob. p4 the die lands on side 4, and
# at each die roll, you keep track of the number of appearances of each side,
# the game stops until either the number of appearances of side 1 equals k,
# OR the number of appearances of side 2 equals k,
# OR the number of appearances of side 3 equals k,
# OR the number of appearances of side 4 equals k, where 1<=k<=n. Try:
# ExpG(10,1/4,1/4,1/4,1/4);
ExpG:=proc(n,p1,p2,p3,p4) local t,k:
[seq(subs(t=1,diff(G(k,k,k,k,p1,p2,p3,p4,t),t)),k=1..n)]
end:



#Seqs1(M,i): Given a list M=[m1,m2,...,mk] and a location 1<=i<=k,
# generates all sequences s of length nops(M) where s[i]=M[i] and s[j]<M[j] j<>i
Seqs1:=proc(M,i) local m,i1,T1,T,t1,j:
option remember:
m:=nops(M):
if m=1 then
 if i=1 then
  return {[M[1]]}:
 else
  return {seq([i1],i1=0..M[1]-1)}:
 fi:
fi:
T:={}:
if m = i then
 T1:=Seqs1([op(1..m-1, M)], i):
 T:={seq([op(t1), M[m]], t1 in T1)}:
else
 T1:=Seqs1([op(1..m-1, M)], i):
 for j from 0 to M[m]-1 do
   T:=T union {seq([op(t1), j], t1 in T1)}:
 od:
fi:
T:
end:

#Seqs(M): the union of all the sequences from Seqs1(M,i) for 1<=i<=nops(M)
Seqs:=proc(M) local i,m:
m:=nops(M):
[seq(Seqs1(M,i),i=1..m)]:
end:


#Fgen(M,P,t): inputs a list M and a list of probabilities P, each of length k,
# and a variable t, outputs the probability generating function of the duration of the game in which
# with probability P[i] the die lands on side i for 1<=i<=k and the game stops until
# the number of appearances of side 1 equals M[1]
# OR the number of appearances of side 2 equals M[2],
# OR the number of appearances of side 3 equals M[3],...,
# OR the number of appearances of side k equals M[k]. Try:
# Fgen([3,3,3],[1/3,1/3,1/3],t);
Fgen:=proc(M,P,t) local S,co,i,k,m,s,j,s1,p:
S:=Seqs(M):
m:=nops(M):
k:=nops(S):
p:=nops(P):
if p<>m then
 print(`Both inputs must be of the same length`):
 return(FAIL):
fi:
 
co:=0:
for i from 1 to k do
 for s in S[i] do
  co:=co+(add(s1,s1 in s)-1)!/mul(s[j]!,j=1..i-1)/(s[i]-1)!/mul(s[j]!,j=i+1..nops(s))*mul(P[j]^s[j],j=1..m)*t^add(s1,s1 in s):
 od:
od:
co:
end:

#ExpFgen(n,P): outputs the first n terms of the expected duration of the game
# in which with probability P[i] the die lands on side i for 1<=i<=k for each k<=n
# and at each die roll, you keep track of the number of appearances of each side,
# the game stops until either the number of appearances of side 1 equals k,
# OR the number of appearances of side 2 equals k,
# OR the number of appearances of side 3 equals k,
# OR the number of appearances of side 4 equals k, where 1<=k<=n. Try:
ExpFgen:=proc(n,P) local t,k,m:
m:=nops(P):
[seq(subs(t=1,diff(Fgen([k$m],P,t),t)),k=1..n)]
end:




#####end of (faster) special cases




#####start Findrec 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:
end:


 
#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:
end:
 

with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:
end:

#####end Findrec