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


print(`First Written: Oct. 2025: tested for Maple 20  `):
print():
print(`This is TFP.txt, a Maple package, accompanying the article:`):
print(`Multi-Gender Talmudic Family Planning`):
print(`by Lucy Martinez and Doron Zeilberger `):


print(`This Maple package, as well as the article, are `):
print(`available from:`):
print(` http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/tfp.html `):
print(``):
print(`Please report bugs to DoronZeil@gmail.com `):

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

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

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

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


print(`For a list of the Pre-computed recurrence operators procedures type ezraOpe(), for help with`):
print(`a specific procedure, type ezra(procedure_name)`):
print(`-----------------------------`):
print(``):
print(``):



ezraOpe:=proc()
if args=NULL then
print(`The opeartor procedures are: `):
print(` OpeAve2, Ope4same, OpeFave, OPEAveF `):
else
ezra(args):
fi:
end:

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

ezra1:=proc()
if args=NULL then
print(`The supporting functions are: CheckP, Chop, comp, easyChop, easycomp, generalcomp, generalChop`):
print(`f1, F1, f2, F2`):
else
ezra(args):
fi:
end:



ezra:=proc()

if args=NULL then
print(`The main procedures are: CheckEasyGenderRatioConjecture, CheckGenderRatioConjecture, easyAveD, EasyAveSeq, easyGFB, easyGFBd, easyGFBt `):
print(`f, F, Fave, GFB, GFBC, GFBCt, GFBCgeneral, GFBgeneral, GFBCgeneralT,  `):
print(`Stat, StatEasyFS, StatFS `):


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]=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. 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 size of G,`):
print(`and a list G with the number of genders, outputs the part with exponents`):
print(`that satisfy reaching the goal of having each of the genders 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]=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]=EasyAveSeq then 
print(`EasyAveSeq(P,K), You have k:=nops(P), genders, and each of them has prob. P[i], outputs the sequence of average family size if your goal is to have at least i of some gender, for i from 1 to K. Try:`):
print(`EasyAveSeq([1/3,2/3],50);`):

elif nargs=1 and args[1]=easyAveD then 
print(`easyAveD(G,P): A direct way to compute the expected number of babies until achieving ONE of the goals in G with prob. dist. P. Try:`):
print(`easyAveD([5,5],[1/2,1/2]);`):

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 size of G,`):
print(`and a list G with the number of genders, outputs the part with exponents`):
print(`that satisfy reaching the goal of having one of the genders in 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 each gender,`):
print(`a list of probabilities P (for each gender), a variable x and`):
print(`outputs the prob. generating function in x[1],...,x[k] where k is the size of G`):
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 babies. 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]=f then 
print(`f(B,G,p,q,x1,x2): The weight-enumerator over all scenarios with B boys and G girls until EITHER reaching B boys or G girls. . Try:`):
print(`f(4,6,1/3,2/3,x1,x2); `):

elif nargs=1 and args[1]=F then 
print(`F(B,G,p,q,x1,x2): The weight-enumerator over all scenarios with B boys and G girls. Try:`):
print(`F(4,6,1/3,2/3,x1,x2); `):


elif nargs=1 and args[1]=Fave then 
print(`Fave(B,G,p): The expected family size for reaching at least B boys and at least G girls woth prob(B)=p and prob(G)=1-p. Try:`):
print(`Fave(4,6,1/3); `):

elif nargs=1 and args[1]=f1 then 
print(`f1(B,G,p,q,x1,x2): The weight-enumerator over all scenarios with B boys and G girls until you either reached B boys or G girls, and reached G first. Try:`):
print(`f1(4,6,1/3,2/3,x1,x2); `):

elif nargs=1 and args[1]=F1 then 
print(`F1(B,G,p,q,x1,x2): The weight-enumerator over all scenarios with B boys and G girls where the G goal was achieved first. Try:`):
print(`F1(4,6,1/3,2/3,x1,x2); `):


elif nargs=1 and args[1]=f2 then 
print(`f2(B,G,p,q,x1,x2): The weight-enumerator over all scenarios with B boys and G girls until you either reached B boys or G girls, and reached B first. Try:`):
print(`f1(4,6,1/3,2/3,x1,x2); `):

elif nargs=1 and args[1]=F2 then 
print(`F2(B,G,p,q,x1,x2): The weight-enumerator over all scenarios with B boys and G girls where the B goal was achieved first. Try:`):
print(`F2(4,6,1/3,2/3,x1,x2); `):

 
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 size of G,`):
print(`and a list G with the number of each gender, outputs the part with exponents`):
print(`that satisfy reaching the goal of having EXACTLY g of the genders in 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 genders,`):
print(`a list of probabilities P (for each gender), a variable x, and an integer K`):
print(`outputs the prob. generating function in x[1],...,x[k] where k is the size of G`):
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] first) or a[2]=g[2] (if you reached G[2] first), etc.`):
print(`assuming that you finished in <=K steps and also`):
print(`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 taking a total of <=K children. 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 children where your goal is to have EXACTLY g of the genders in 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 babies. 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): 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 babies. Try:`):
print(`GFBCt([5,3,2],[1/2,1/8,3/8],t,10); `):


elif nargs=1 and args[1]=GFBgeneral then 
print(`GFBgeneral(G,P,x,K,g): inputs a list G corresponding to the number of each gender,`):
print(`a list of probabilities P (for each gender), 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 size of G`):
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 of the genders in the list G, the prob.`):
print(`either a[1]=g[1] (if you reached G[1] first) or a[2]=g[2] (if you reached G[2] first), etc.`):
print(`assuming that you finished in <=K steps 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]=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 midSimuG(L,P) 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]=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): outputs the operator that satisfies L:=EasyAveSeq([p,1-p],50); where `):
print(`where L is the sequence of average family size if your goal is to have at least i of some gender, where 1<=i<=50. Try:`):
print(`OpeAve2(n,N,1/3); `):


elif nargs=1 and args[1]=OPEAveF then 
print(`OPEAveF(n,m,N,M,p): the paris of operators [opeN,opeM] such that`):
print(`AveF(n,m)=OpeN Ave(n,m) and AveF(n,m)=OpeM Ave(n,m) . Try:`):
print(`OPEAveF(n,m,N,M,p); `):


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


elif nargs=1 and args[1]=OpeFave then 
print(`OpeFave(n,N,p):outputs the operator that annihilates Fave(n,n,p). Try:`):
print(`OpeFave(n,N,1/3);`):

elif nargs=1 and args[1]=Simu then 
print(`Simu(p,b,g): Inputs a rational number 0<=p<=1, and pos. integers b and g`):
print(`and with probability p gives a boy and with probility 1-p gives a girl`):
print(`and keeps going until you have reached your goal of having both (at least)`):
print(`b boys and (at least) g girls. The output is a pair [m,n] where`):
print(`m is the number of boys you have and n is the number of girls you have. 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);`):


else
print(`There is no ezra for`,args):
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,b,g): inputs a rational number 0<=p<=1, and pos. integers b and g
# and with probability p gives a boy and with probility 1-p gives a girl
# and keeps going until you have reached your goal of having both (at least)
# b boys and (at least) g girls. The output is a pair [m,n] where
# m is the number of boys you have and n is the number of girls you have
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,p,X,X1:
k:=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 k<>nops(L) then
 return(FAIL):
fi:
X:=[0$k]:

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 genders,
# a list of probabilities P (for each gender), a variable x, and an integer K
# outputs the prob. generating function in x[1],...,x[k] where k is the size of G
# 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] first) or a[2]=g[2] (if you reached G[2] first), etc.
# assuming that you finished in <=K steps 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 size of G,
# and a list G with the number of genders, outputs the part with exponents
# that satisfy reaching the goal of having each of the genders 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 taking a total of <=K children
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 babies
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:


#F1(B,G,p,q,x1,x2): The weigt-enumerator over all scenarios with B boys and G girls where the G goal was achieved first. Try:
#F1(4,6,1/3,2/3,x1,x2);
F1:=proc(B,G,p,q,x1,x2) local b1,g1:
sum(sum(binomial(G+b1-1,G-1)*binomial(B-1-b1+g1-G,B-1-b1)*(p*x1)^B*(q*x2)^g1,g1=G..infinity),b1=0..B-1):
end:


#F1ave(B,G,p,q): The expected family size over all scenarios with B boys and G girls where the G goal was achieved first. Try:
#F1ave(4,6,1/3,2/3);
F1ave:=proc(B,G,p,q,x1,x2) local b1,g1:
sum(sum(binomial(G+b1-1,G-1)*binomial(B-1-b1+g1-G,B-1-b1)*(p)^B*(q)^g1*(B+g1),g1=G..infinity),b1=0..B-1):
end:


#F2(B,G,p,q,x1,x2): The weigt-enumerator over all scenarios with B boys and G girls where the B goal was achieved first. Try:
#F2(4,6,1/3,2/3,x1,x2);
F2:=proc(B,G,p,q,x1,x2) local b1,g1:
sum(sum(binomial(B+g1-1,B-1)*binomial(G-1-g1+b1-B,G-1-g1)*(p*x1)^b1*(q*x2)^G,b1=B..infinity),g1=0..G-1):
end:


#F2ave(B,G,p,q): The expected family zize over all scenarios with B boys and G girls where the B goal was achieved first. Try:
#F2ave(4,6,1/3,2/3);
F2ave:=proc(B,G,p,q) local b1,g1:
sum(sum(binomial(B+g1-1,B-1)*binomial(G-1-g1+b1-B,G-1-g1)*(p)^b1*(q)^G*(b1+G),b1=B..infinity),g1=0..G-1):
end:

#F(B,G,p,q,x1,x2): The weigt-enumerator over all scenarios with B boys and G girls. Try:
#F(4,6,1/3,2/3,x1,x2);
F:=proc(B,G,p,q,x1,x2):
F1(B,G,p,q,x1,x2)+F2(B,G,p,q,x1,x2):
end:


#Fave(B,G,p): The expected family size over all scenarios with B boys and G girls. Try:
#Fave(4,6,1/3);
Fave:=proc(B,G,p):
F1ave(B,G,p,1-p)+F2ave(B,G,p,1-p):
end:



#f1(B,G,p,q,x1,x2): The weigt-enumerator over all scenarios until you either reached B boys or reached G girls and reached B first. Try:
#f1(4,6,1/3,2/3,x1,x2);
f1:=proc(B,G,p,q,x1,x2) local b1:
sum(binomial(G+b1-1,G-1)*(p*x1)^b1*(q*x2)^G,b1=0..B-1):
end:


#f2(B,G,p,q,x1,x2): The weigt-enumerator over all scenarios until you either reached B boys or reached G girls and reached G first. Try:
#f2(4,6,1/3,2/3,x1,x2);
f2:=proc(B,G,p,q,x1,x2) local g1:
sum(binomial(B-1+g1,B-1)*(p*x1)^B*(q*x2)^g1,g1=0..G-1):
end:



#f(B,G,p,q,x1,x2): The weigt-enumerator over all scenarios with B boys and G girls until you either reached B first or G first. Try:
#f(4,6,1/3,2/3,x1,x2);
f:=proc(B,G,p,q,x1,x2):
f1(B,G,p,q,x1,x2)+f2(B,G,p,q,x1,x2):
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,p,X,X1:
k:=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 k<>nops(L) then
 return(FAIL):
fi:
X:=[0$k]:

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 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:=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 each genders,
# a list of probabilities P (for each gender), a variable x,
# outputs the prob. generating function in x[1],...,x[k] where k is the size of G
# 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 size of G,
# and a list G with the number of genders, outputs the part with exponents
# that satisfy reaching the goal of having one of the genders in 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 babies
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:


#EasyAveSeqOld(P,K), You have k:=nops(P), genders, and each of them has prob. P[i], outputs the sequence of average family size if your goal is to have at least i of some gender, for i from 1 to K. Try:
#EasyAveSeqOld([1/3,2/3],50);
EasyAveSeqOld:=proc(P,K) local i,L,k,t:
k:=nops(P):
L:=[seq(easyGFBt([i$k],P,t),i=1..K)]:
[seq(subs(t=1,diff(L[i],t)),i=1..K)]:
end:


#EasyAveSeq(P,K), You have k:=nops(P), genders, and each of them has prob. P[i], outputs the sequence of average family size if your goal is to have at least i of some gender, for i from 1 to K. Try:
#EasyAveSeq([1/3,2/3],50);
EasyAveSeq:=proc(P,K) local k,i:
k:=nops(P):
[seq(easyAveD([i$k],P),i=1..K)]:
end:



#OpeAve2(n,N,p): outputs the operator that satisfies L:=EasyAveSeq([p,1-p],50);
# where L is the sequence of average family size if your goal is to have at least i of some gender, where 1<=i<=50.
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:

#OpeFaveSame(n,N):outputs the operator that annihilates Fave(1/2,1/2,n,n). Try:
#OpeFaveSame(n,N);
OpeFaveSame:=proc(n,N):1/2*(2*n+1)/n-1/2*(4*n+5)/(n+1)*N+N^2:end:

#OpeFave(n,N,p):outputs the operator that annihilates Fave(p,1-p,n,n). Try:
#OpeFave(n,N,p);
OpeFave:=proc(n,N,p):
-2*(-1+p)*(2*n+1)*p/n+(4*n*p^2-4*n*p+2*p^2-n-2*p-2)/(n+1)*N+N^2: end:



#Ope4same(n,N): the linear recurrence operator for the sequence
#of expected number of babies with 4 genders each with prob. 1/4
#you finish as soon as one of genders has n babies. Try:
#Ope4same(n,N):

Ope4same:=proc(n,N):
-1/576*(4*n+1)*(3*n+2)*(2*n+3)*(2*n+1)*(3*n+1)*(4*n+3)*(5514198858742237*n^2+
22377255642631594*n+21709635674546661)/(n+2)/(28152901853008477*n^2+
81763683068059260*n+53615706600102264)/(n+3)^2/(n+4)^3+1/576*(2*n+3)*(
14460392818937498400*n^7+182606360698176085584*n^6+1052157967879148282152*n^5+
3484773784146230875959*n^4+6979199009322021232333*n^3+8314830354132186721782*n^
2+5412609663886380033970*n+1480437544492104202620)/(n+2)/(28152901853008477*n^2
+81763683068059260*n+53615706600102264)/(n+3)^2/(n+4)^3*N-1/576*(
83921357125144702080*n^8+1443634929829166314752*n^7+11054199574946609460984*n^6
+48676841768444583401754*n^5+133788726394565359658831*n^4+
233751077270618823645175*n^3+252381873231581655737525*n^2+
153212526518476108150989*n+39775428194986725479910)/(n+2)/(28152901853008477*n^
2+81763683068059260*n+53615706600102264)/(n+3)^2/(n+4)^3*N^2+1/576*(
110001142974539410560*n^7+1915279238626411246848*n^6+14304392985937467303200*n^
5+58976182566812757393778*n^4+144026256606834288708225*n^3+
206839081177665128664287*n^2+160221387224420715049270*n+50931907675076176759332
)/(28152901853008477*n^2+81763683068059260*n+53615706600102264)/(n+3)^2/(n+4)^3
*N^3-1/288*(34020232205983529760*n^5+450870038338456910016*n^4+
2327471141692767026974*n^3+5760259258259629592377*n^2+6685431166857223519260*n+
2835874645996633400088)/(28152901853008477*n^2+81763683068059260*n+
53615706600102264)/(n+4)^3*N^4+N^5:
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:


#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
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:



#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,p,X,X1:
k:=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 k<>nops(L) then
 return(FAIL):
fi:
X:=[0$k]:

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 midSimuG(L,P) and divides by K
ManygeneralSimuG:=proc(L,P,K,x,g) local ele,i,k,poly,M,i1:
#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[i1]^M[i1],i1=1..k):
od:
poly/K:
end:


#GFBgeneral(G,P,x,K,g): inputs a list G corresponding to the number of each gender,
# a list of probabilities P (for each gender), 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 size of G
# 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) genders, the prob.
# either a[1]=g[1] (if you reached G[1] first) or a[2]=g[2] (if you reached G[2] first), etc.
# assuming that you finished in <=K steps 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):
 f1:=generalChop(F,x,G,g):
 AlreadyDone:=expand(AlreadyDone+f1):
 F:=F-f1:
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 size of G,
# and a list G with the number of each gender, outputs the part with exponents
# that satisfy reaching the goal of having exactly g of the genders in the list G
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 children
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 babies
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:

####end of general


#start direct

#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:


#easyAveD(G,P): A direct way to compute the expected number of babies until achieving ONE of the goals in G with prob. dist. O. Try:
#easyAveD([5,5],[1/2,1/2]);

easyAveD:=proc(G,P) 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])^G[i]*
add(
(G[i]+convert(v,`+`))*
mul((P[i1])^v[i1],i1=1..i-1)*mul((P[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:



#end direct



#OPEAveF(n,m,N,M,p): the paris of operators [opeN,opeM] such that
#AveF(n,m)=OpeN Ave(n,m) and AveF(n,m)=OpeM Ave(n,m) . Try:
#OPEAveF(n,m,N,M,p):
OPEAveF:=proc(n,m,N,M,p):
[-(m*p+n*p-3*m-n-2*p+4)/(m-1)/M+(2*m*p+2*n*p-3*m-2*n-4*p+5)/(m-1)/M^2-(-1+p)*(m-2+n)/(m-1)/M^3, (m*p+n*p+2*n-2*p-2)/(n-1)/N-(2*m*p+2*n*p+n-4*p-1)/(n-1)/N^2+p*(m-2+n)/(n-1)/N^3]
end:



FaveF:=proc(n1,m1,p) local ope,n,m,N,M,i1,L:
option remember:
ope:=
[-(m*p+n*p-3*m-n-2*p+4)/(m-1)/M+(2*m*p+2*n*p-3*m-2*n-4*p+5)/(m-1)/M^2-(-1+p)*(m-2+n)/(m-1)/M^3, (m*p+n*p+2*n-2*p-2)/(n-1)/N-(2*m*p+2*n*p+n-4*p-1)/(n-1)/N^2+p*(m-2+n)/(n-1)/N^3]:


L:=
[[-p+1/p-1/(1-p)*p^2+2/(1-p)*p, 1/p-1/(1-p)*p^2+3/(1-p)*p, -p^2+2*p+1/p-1/(1-p)*p^2+4/(1-p)*p], [2/p-2*p^2-2/(1-p)*p^3+3/(1-p)*p^2, 2/p-2/(1-p)*p^3+4/(1-p)*p^2, -3*p^3+5*p^2+2/p-2/(1-p)*p^3+5/(1-p)*p^2], [3/p-3*p^3-3/(1-p)*p^4+4/(1-p)*p^3, 3/p-3/(1-p)*p^4+5/(1-p)*p^3, 9*p^
3+3/p-6*p^4-3/(1-p)*p^4+6/(1-p)*p^3]]:

if n1<=3 and  m1<=3 then
 RETURN(L[n1][m1]):
 fi:

if n1<=3 and m1>=3  then
RETURN(expand(add(subs({n=n1,m=m1},coeff(ope[1],M,-i1))*FaveF(n1,m1-i1,p),i1=1..3))):
fi:


if n1>3  then
RETURN(expand(add(subs({n=n1,m=m1},coeff(ope[2],N,-i1))*FaveF(n1-i1,m1,p),i1=1..3))):
fi:

FAIL:
end: