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

print(`First Written :Nov. 2025: tested for Maple 20  `):
print():
print(`This is CoinToss.txt, a Maple package, accompanying the article:`):
print(` How many coin tosses does it take until you get r Heads OR s Tails?`):
print(`by Svante Janson, 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/cointoss.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(`For a list of the OLD procedures type ezraOld(), 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 story  procedures type ezraStory(), 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(`For a list of the Fast procedures type ezraFast(), for help with`):
print(`a specific procedure, type ezra(procedure_name)`):
print(`-----------------------------`):
print(``):
print(``):



ezraOld:=proc()
if args=NULL then
print(`The Old procedures are: `):
print(`CheckEasyGenderRatioConjecture, CheckGenderRatioConjecture, easyAveD, EasyAveSeq, easyGFB, easyGFBd, easyGFBt `):
print(`GFB, GFBC, GFBCt, GFBCgeneral, GFBgeneral, GFBCgeneralT,`):
else
ezra(args):
fi:
end:


ezraStory:=proc()
if args=NULL then
print(`The story procedures are:  FairMomStory, TFPstoryGFbi, TFPstoryE, TFPstoryGFd `):
print(`  `):
else
ezra(args):
fi:
end:

ezraFast:=proc()
if args=NULL then
print(`The FAST procedures (using recurrences via WZ theory) are: `):
print(` FairFacMomF, FairMomF, FairMomCF1, ftDfast, faveFd, FtDfast, FaveF, FaveFd, Ffast, faveF, ffast `):
else
ezra(args):
fi:
end:

ezraOpe:=proc()
if args=NULL then
print(`The opeartor procedures are: `):
print(` OpefAve, Ope4same, OpeFave, OPEAveF,  OpeFt, OPEx1x2 `):
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: AbsNorMoms, AbsNorMomsN, CheckP, Chop, comp, easyChop, easycomp, generalcomp, generalChop`):
print(`CheckEB, f1, F1, f2, F2, FairFacMom, plotdist, rf`):
else
ezra(args):
fi:
end:



ezra:=proc()

if args=NULL then
print(`The main procedures are: `):
print(`AveAndMoms2S, f, fave, F, Fave,  FairMom, FairMomsC, Stat, StatEasyFS, StatFS `):


elif nargs=1 and args[1]=AbsNorMoms then
print(`AbsNorMoms(R)`):
print(`The mean, variance and the Scaled centered moments of |N(0,1)|. (pdf exp(-x^2/2)/sqrt(Pi/2) supported in [0,infinity].) Try: `):
print(`AbsNorMoms(6);`):

elif nargs=1 and args[1]=AbsNorMomsN then
print(`AbsNorMomsN(R)`):
print(`The mean, variance and the Scaled centered moments of -|N(0,1)|.`):
print(`(pdf exp(-x^2/2)/sqrt(Pi/2) supported in [-infinity,0].). Try: `):
print(`AbsNorMomsN(6);`):

elif nargs=1 and args[1]=AveAndMoms2 then 
print(`AveAndMoms2(f,x1,x2,N): Given a specific`):
print(`probability generating function`):
print(`for a bivariate distribution`):
print(`f(x1,x2) (a polynomial, rational function and possibly beyond)`):
print(`returns the averages of the first and second r.v. followed`):
print(`by a list-of-lists L`):
print(`such that L[r+1][s+1] is the (r,s)-mixed moment about the mean`):
print(`For example, try:`):
print(`AveAndMoms2(((1+x1+x2)/3)^100,x1,x2,4);`):


elif nargs=1 and args[1]=AveAndMoms2S then 
print(`AveAndMoms2S(f,x1,x2,N): Given a specific`):
print(`probability generating function`):
print(`for a bivariate distribution`):
print(`f(x1,x2) (a polynomial, rational function and possibly beyond)`):
print(`returns the averages of the first and second r.v. followed`):
print(`by a list-of-lists L`):
print(`such that L[r+1][s+1] is the SCALED (r,s)-mixed moment about the mean`):
print(`For example, try:`):
print(`AveAndMoms2S(((1+x1+x2)/3)^100,x1,x2,4);`):


elif nargs=1 and args[1]=CheckEB then 
print(`CheckEB(a,b,R): Checks the explicit formulas for  Fave(a*n,b*n,a/(a+b)) and fave(a*n,b*n,a/(a+b)) for n from 1 to R. Try:`):
print(`CheckEB(3,11,20);`):

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 Heads and G Tails until EITHER reaching B Heads or G girls. . Try:`):
print(`f(4,6,1/3,2/3,x1,x2); `):

elif nargs=1 and args[1]=FairFacMom then 
print(`FairFacMom(n,r): The r-th Factorial moment of the number of coin-tosses of a fair coin until reaching n Heads or n Tails. Try:`):
print(`FairFacMom(n,2);`):

elif nargs=1 and args[1]=FairMom then 
print(`FairMom(n,r): The r-th moment of the number of coin-tosses of a fair coin until reaching n Heads or n Tails. Try:`):
print(`FairMom(n,2);`):

elif nargs=1 and args[1]=FairMomsC then 
print(`FairMomsC(n,R): The expectation, variance, and central moments for reaching n Heads for n Tails with a fair coin up to the R-th moment.`):
print(`It also returns the limits of the scaled moments and the floating point of the latter.Try:`):
print(`FairMomsC(n,4);`):


elif nargs=1 and args[1]=FairMomCF1 then 
print(`FairMomsCF1(n,R,C): The expectation, variance, and central moments for reaching n Heads for n Tails with a fair coin up to the R-th moment.`):
print(`in terms of C:=binomial(2*n,n+1)*(n+1)/4^n:`):
print(`It also returns the limits of the scaled moments and the floating point of the latter.Try:`):
print(`FairMomsCF1(n,4,C);`):


elif nargs=1 and args[1]=FairMomStory then 
print(`FairMomStory(n,R,S): An article about the time it takes to reach n Heads or n Tails with a fair coin.`):
print(`It gives the explicit expressions for the moments, and the scaled limits from 1 to R`):
print(`and for those from R+1 to S just checks that the limit of the scaled central moments are`):
print(`the same as those of -|N(0,1)| Try:`):
print(`FairMomStory(n,4,20);`):


elif nargs=1 and args[1]=fave then 
print(`fave(B,G,p): The expected family size for reaching at least B Heads OR at least G girls (whaever comes first). With prob(B)=p and prob(G)=1-p. Try:`):
print(`fave(4,6,1/3); `):


elif nargs=1 and args[1]=faveF then 
print(`faveF(B,G,p): same as fave(B,G,p) but done faster, using the recurrences.`):
print(`faveF(4,6,1/3); `):


elif nargs=1 and args[1]=faveFd then 
print(`faveFd(n,p): A fast way to compute fave(n,n,p), using the 2nd order recurrences. Try:`):
print(`faveFd(10,p);`):

elif nargs=1 and args[1]=ffast then 
print(`ffast(B,G,p,x1,x2): The weight-enumerator over all scenarios with B Heads and G girls until EITHER reaching B Heads or G girls. Done FAST, using the recurrences . Same as f(B,G,p,1-p,x1,x2), but faster.Try:`):
print(`ffast(4,6,1/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 Heads and G girls. Try:`):
print(`F(4,6,1/3,2/3,x1,x2); `):


elif nargs=1 and args[1]=Ffast then 
print(`Ffast(B,G,p,x1,x2): The weight-enumerator over all scenarios with B Heads and G girls with Probability of B being p (and hence prob. of girl being 1-p). DONE FAST (via the recurrences) Try:`):
print(`Same as F(B,G,p,1-p,x1,x2)`):
print(`Ffast(4,6,1/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 Heads AND at least G girls with prob(B)=p and prob(G)=1-p.done directly (so slow), Try:`):
print(`Fave(4,6,1/3); `):


elif nargs=1 and args[1]=FaveF then 
print(`FaveF(n,m,p): A fast way to compute Fave(n,m,p), using the pure 3rd order recurrences in each direction. Try:`):
print(`FaveF(10,10,p);`):

elif nargs=1 and args[1]=FaveFd then 
print(`FaveFd(n,p): A fast way to compute Fave(n,n,p), using the 2nd order recurrences. Try:`):
print(`FaveFd(10,p);`):


elif nargs=1 and args[1]=f1 then 
print(`f1(B,G,p,q,x1,x2): The weight-enumerator over all scenarios with B Heads and G girls until you either reached B Heads 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 Heads 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 Heads and G girls until you either reached B Heads 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 Heads 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]=ftDfast then 
print(`ftDfast(p,t,n): a fast way to compute f(n,n,p,1-p,t,t):  Try:`):
print(`[seq(ftDfast(p,t,i),i=1..10)];`):


elif nargs=1 and args[1]=FtDfast then 
print(`FtDfast(p,t,n): a fast way to compute F(n,n,p,1-p,t,t):  Try:`):
print(`[seq(FtDfast(p,t,i),i=1..10)];`):

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]=OpefAve then 
print(`OpefAve(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(`OpefAve(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]=OpeFt then 
print(`OpeFt(p,t,n,N): The third-order recurrence operator annihilating  BOTH F(n,n,p,1-p,t,t) and f(n,n,p,1-p,t,t) . `):
print(`Try:`):
print(`OpeFt(p,t,n,N);`):

elif nargs=1 and args[1]=OPEx1x2 then 
print(`OPEx1x2(p,x1,x2,n,m,N,M): The pair of pure operators in n and m (with shift operators N and M) annihilating both`):
print(`f(n,m,p,1-p,x1,x2) and F(n,m,p,1-p,x1,x2) that enable fast computations. Try:`):
print(`OPEx1x2(p,x1,x2,n,m,N,M);`):


elif nargs=1 and args[1]=plotDist then 
print(`plotDist(f,x,K): Given a prob. gen. function f(x) that has a `):
print(`Taylor series`):
print(`for a discrete r.v.`):
print(`plots its normalized version (X-mu)/sig between mu-K1*sig and`):
print(`mu+K2*sig`):
print(`For example, try:`):
print(`plotDist((1+x)^40,x,5,5);`):

elif nargs=1 and args[1]=rf then 
print(`rf(a,n): The raising factorial a(a+1)...(a+n-1). Try:`):
print(`rf(1/2,7);`):

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 Heads and (at least) g girls. The output is a pair [m,n] where`):
print(`m is the number of Heads 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);`):


elif nargs=1 and args[1]=TFPstoryGFd then 
print(`TFPstoryGFd(p,t,K1,K2): inputs a numeric or symbolic prob. p a symbol t, and integer K outputs an article about the `):
print(`prob. generating functions for the famiy sizes until reaching n Heads AND n girls, it gives the`):
print(`recurrence (the same) for the case where you quit as soon as you have n babies of one gender and`):
print(`the case where you need at least n babies from both genders.`):
print(`then it compares the statistics of these two kinds of goals`):
print(`TFPstoryGFd(p,t,10,4);`):

elif nargs=1 and args[1]=TFPstoryE then 
print(`TFPstoryE(p,K): inputs a numeric or symbolic prob. p and integer K outputs an article about the quantity`):
print(`Expected family size in Talmudic Family Planning until reaching B Heads AND G girls, it gives the`):
print(`two pure recurrences and a K by K matrix (as a list of lists) such that M[B][G] is the expected`):
print(`family size if you goal is to have B Heads and G girls. Try:`):
print(`TFPstoryE(p,10);`):


elif nargs=1 and args[1]=TFPstoryGFbi then 
print(`TFPstoryGFbi(p,K)): inputs a numeric or symbolic prob. p and integer K outputs an article about the bi-variate`):
print(`prob. generating function for the number of heads and the number of tails until reaching either`):
print(`h Heads OR t TAILS`):
print(`and also when the goal is`):
print(`h Heads AND t TAILS`):
print(`where Pr(H)=p and Pr(T)=1-p`):
print(`It gives the linear recurrences (the SAME) and the respective initial conditions and`):
print(`then does some statistical analysis. Try:`):
print(`TFPstoryGFbi(p,4);`):

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 Heads and (at least) g girls. The output is a pair [m,n] where
# m is the number of Heads 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:


#F1old(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:
#F1old(4,6,1/3,2/3,x1,x2);
F1old:=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:

#F1old(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:
#F1old(4,6,1/3,2/3,x1,x2);
F1old:=proc(B,G,p,q,x1,x2) local g:
sum(binomial(B-1+g,B-1)*(p*x1)^B*(q*x2)^g,g=G..infinity):
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 g:
normal(((p*x1)/(1-q*x2))^B-f2(B,G,p,q,x1,x2)):
end:


#F1aveOld(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:
#F1aveOld(4,6,1/3,2/3);
F1aveOld:=proc(B,G,p,q) local b1,g1:
sum((B+g)*binomial(B-1+g,B-1)*(p)^B*(q)^g,g=G..infinity):
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) local b1,g1:
normal(p^B*B/(1-q)^(B+1)-f2ave(B,G,p,q,x1,x2)):
end:


#F2old(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:
#F2old(4,6,1/3,2/3,x1,x2);
F2old:=proc(B,G,p,q,x1,x2) local b    :
sum(binomial(G-1+b,G-1)*(p*x1)^b*(q*x2)^G,b=B..infinity):
end:

F2:=proc(B,G,p,q,x1,x2) local g:
normal(((q*x2)/(1-p*x1))^G-f1(B,G,p,q,x1,x2)):
end:


#F2aveOld(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:
#F2aveOld(4,6,1/3,2/3);
F2aveOld:=proc(B,G,p,q) local b1,g1:
sum((b+G)*binomial(G-1+b,G-1)*(p)^b*(q)^G,b=B..infinity):
end:

F2ave:=proc(B,G,p,q) local b1,g1:
normal(q^G*G/(1-p)^(G+1)-f1ave(B,G,p,q,x1,x2)):
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):
normal(F1ave(B,G,p,1-p)+F2ave(B,G,p,1-p)):
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:



#f1ave(B,G,p,q): The average over all scenarios until you either reached B boys or reached G girls and reached B first. Try:
#f1ave(4,6,1/3,2/3);
f1ave:=proc(B,G,p,q) local b1:
sum((b1+G)*binomial(G+b1-1,G-1)*(p)^b1*(q)^G,b1=0..B-1):
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:


#f2ave(B,G,p,q): The average over all scenarios until you either reached B boys or reached G girls and reached G first. Try:
#f2ave(4,6,1/3,2/3);
f2ave:=proc(B,G,p,q) local g1:
sum((B+g1)*binomial(B-1+g1,B-1)*(p)^B*(q)^g1,g1=0..G-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:=F/subs(t=1,F):
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:


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:



#OpefAve(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.
OpefAve:=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:=
normal([[-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(normal(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(normal(add(subs({n=n1,m=m1},coeff(ope[2],N,-i1))*FaveF(n1-i1,m1,p),i1=1..3))):
fi:

FAIL:
end:

#FaveFd(n,p): A fast way to compute Fave(n,n,p). Try:
#FaveFd(10,p);
FaveFd:=proc(n,p) option remember:
if n=1 then
-(p^2-p+1)/(-1+p)/p
elif n=2 then
2*(p^4-2*p^3+p-1)/(-1+p)/p:
else
normal(-(4*n*p^2-4*n*p-6*p^2-n+6*p)/(n-1)*FaveFd(n-1,p)+
2*(p-1)*(2*n-3)*p/(n-2)*FaveFd(n-2,p)):
fi:
end:


#faveFd(n,p): A fast way to compute fave(n,n,p). Try:
#faveFd(10,p);
faveFd:=proc(n,p) option remember:
if n=1 then
1:
elif n=2 then
-2*p^2+2*p+2
else
normal(-(4*n*p^2-4*n*p-6*p^2-n+6*p)/(n-1)*faveFd(n-1,p)+
2*(p-1)*(2*n-3)*p/(n-2)*faveFd(n-2,p)):
fi:
end:


#OpeFt(p,t,n,N): The third-order recurrence operator annihilating  F(n,n,p,1-p,t,t);
#Try:
#OpeFt(p,t,n,N);
OpeFt:=proc(p,t,n,N)
2*p^2*t^4*(-1+p)^2*(2*n+1)*(4*n*p^2*t^2-4*n*p*t^2+6*p^2*t^2-6*p*t^2+2*n*t-n+3*t
-2)/(n+2)/(4*n*p^2*t^2-4*n*p*t^2+2*p^2*t^2-2*p*t^2+2*n*t-n+t-1)/(p*t-1)/(p*t-t+
1)-p*t^2*(-1+p)*(32*n^2*p^4*t^4-64*n^2*p^3*t^4+64*n*p^4*t^4+48*n^2*p^2*t^4-128*
n*p^3*t^4+24*p^4*t^4-16*n^2*p*t^4+96*n*p^2*t^4-48*p^3*t^4-12*n^2*p^2*t^2-32*n*p
*t^4+36*p^2*t^4+12*n^2*p*t^2+8*n^2*t^3-30*n*p^2*t^2-12*p*t^4-12*n^2*t^2+30*n*p*
t^2+16*n*t^3-12*p^2*t^2+2*n^2*t-26*n*t^2+12*p*t^2+6*t^3+n^2+5*n*t-10*t^2+3*n+2*
t+2)/(n+2)/(4*n*p^2*t^2-4*n*p*t^2+2*p^2*t^2-2*p*t^2+2*n*t-n+t-1)/(p*t-1)/(p*t-t
+1)*N+t*(16*n^2*p^6*t^5-48*n^2*p^5*t^5+32*n*p^6*t^5+48*n^2*p^4*t^5-96*n*p^5*t^5
+12*p^6*t^5+24*n^2*p^4*t^4-16*n^2*p^3*t^5+96*n*p^4*t^5-36*p^5*t^5-28*n^2*p^4*t^
3-48*n^2*p^3*t^4+48*n*p^4*t^4-32*n*p^3*t^5+36*p^4*t^5+56*n^2*p^3*t^3+24*n^2*p^2
*t^4-62*n*p^4*t^3-96*n*p^3*t^4+18*p^4*t^4-12*p^3*t^5-24*n^2*p^2*t^3+124*n*p^3*t
^3+48*n*p^2*t^4-24*p^4*t^3-36*p^3*t^4-12*n^2*p^2*t^2-4*n^2*p*t^3-56*n*p^2*t^3+
48*p^3*t^3+18*p^2*t^4+6*n^2*p^2*t+12*n^2*p*t^2-26*n*p^2*t^2-6*n*p*t^3-22*p^2*t^
3-6*n^2*p*t-2*n^2*t^2+16*n*p^2*t+26*n*p*t^2-10*p^2*t^2-2*p*t^3+3*n^2*t-16*n*p*t
-5*n*t^2+8*p^2*t+10*p*t^2-n^2+8*n*t-8*p*t-2*t^2-3*n+4*t-2)/(n+2)/(4*n*p^2*t^2-4
*n*p*t^2+2*p^2*t^2-2*p*t^2+2*n*t-n+t-1)/(p*t-1)/(p*t-t+1)*N^2+N^3:
end:





#ftDfast(p,t,n): a fast way to compute f(n,n,p,1-p,t,t):  Try:
#[seq(ftDfast(p,t,i),i=1..10)];
ftDfast:=proc(p,t,n1) local L1,A,i1,n:
option remember:
L1:=[t, -2*p^2*t^3+2*p^2*t^2+2*p*t^3-2*p*t^2+t^2, 6*p^4*t^5-6*p^4*t^4-12*p^3*t^5+12*p^3*t^4+6*p^2*t^5-9*p^2*t^4+3*p^2*t^3+3*p*t^4-3*p*t^3+t^3]:

A:=
[-t*(16*n^2*p^6*t^5-48*n^2*p^5*t^5-64*n*p^6*t^5+48*n^2*p^4*t^5+192*n*p^5*t^5+60*p^6*t^5+24*n^2*p^4*t^4-16*n^2*p^3*t^5-192*n*p^4*t^5-180*p^5*t^5-28*n^2*p^4*t^3-48*n^2*p^3*t^4-96*n*p^4*t^4+64*n*p^3*t^5+180*p^4*t^5+56*n^2*p^3*t^3+24*n^2*p^2*t^4+106*n*p^4*t^3+192*n*p^3*t^4+90*
p^4*t^4-60*p^3*t^5-24*n^2*p^2*t^3-212*n*p^3*t^3-96*n*p^2*t^4-90*p^4*t^3-180*p^3*t^4-12*n^2*p^2*t^2-4*n^2*p*t^3+88*n*p^2*t^3+180*p^3*t^3+90*p^2*t^4+6*n^2*p^2*t+12*n^2*p*t^2+46*n*p^2*t^2+18*n*p*t^3-70*p^2*t^3-6*n^2*p*t-2*n^2*t^2-20*n*p^2*t-46*n*p*t^2-40*p^2*t^2-20*p*t^3+3*n^
2*t+20*n*p*t+7*n*t^2+14*p^2*t+40*p*t^2-n^2-10*n*t-14*p*t-5*t^2+3*n+7*t-2)/(n-1)/(4*n*p^2*t^2-4*n*p*t^2-10*p^2*t^2+10*p*t^2+2*n*t-n-5*t+2)/(p*t-1)/(p*t-t+1), p*t^2*(p-1)*(32*n^2*p^4*t^4-64*n^2*p^3*t^4-128*n*p^4*t^4+48*n^2*p^2*t^4+256*n*p^3*t^4+120*p^4*t^4-16*n^2*p*t^4-192*n
*p^2*t^4-240*p^3*t^4-12*n^2*p^2*t^2+64*n*p*t^4+180*p^2*t^4+12*n^2*p*t^2+8*n^2*t^3+42*n*p^2*t^2-60*p*t^4-12*n^2*t^2-42*n*p*t^2-32*n*t^3-30*p^2*t^2+2*n^2*t+46*n*t^2+30*p*t^2+30*t^3+n^2-7*n*t-40*t^2-3*n+5*t+2)/(n-1)/(4*n*p^2*t^2-4*n*p*t^2-10*p^2*t^2+10*p*t^2+2*n*t-n-5*t+2)/(p
*t-1)/(p*t-t+1), -2*p^2*t^4*(p-1)^2*(2*n-5)*(4*n*p^2*t^2-4*n*p*t^2-6*p^2*t^2+6*p*t^2+2*n*t-n-3*t+1)/(n-1)/(4*n*p^2*t^2-4*n*p*t^2-10*p^2*t^2+10*p*t^2+2*n*t-n-5*t+2)/(p*t-1)/(p*t-t+1)]:

if not (type(n1,integer) and n1>0) then
RETURN(FAIL):
fi:

if n1>=1 and n1<=3 then
L1[n1]:
else
normal(add(ftDfast(p,t,n1-i1)*subs(n=n1,A[i1]),i1=1..3)):
fi:

end:




#FtDfast(p,t,n): a fast way to compute F(n,n,p,1-p,t,t):  Try:
#[seq(ftDfast(p,t,i),i=1..10)];
FtDfast:=proc(p,t,n1) local L1,A,i1,n:
option remember:


L1:=
[-p*t^2*(-1+p)*(t-2)/(p*t-t+1)/(p*t-1), p^2*t^4*(-1+p)^2*(2*p^2*t^3-2*p^2*t^2-2*p*t^3+2*p*t^2+3*t^2-8*t+6)/(p*t-t+1)^2/(p*t-1)^2,
-p^3*t^6*(-1+p)^3*(6*p^4*t^5-6*p^4*t^4-12*p^3*t^5+12*p^3*t^4+6*p^2*t^5+9*p^2*t^4-33*p^2*t^3-15*p*t^4+18*p^2*t^2+33*p*t^3-18*p*t^2+10*t^3-36*t^2+45*t-20)/(p*t-t+1)^3/(p*t-1)^3]:

A:=
[-t*(16*n^2*p^6*t^5-48*n^2*p^5*t^5-64*n*p^6*t^5+48*n^2*p^4*t^5+192*n*p^5*t^5+60*p^6*t^5+24*n^2*p^4*t^4-16*n^2*p^3*t^5-192*n*p^4*t^5-180*p^5*t^5-28*n^2*p^4*t^3-48*n^2*p^3*t^4-96*n*p^4*t^4+64*n*p^3*t^5+180*p^4*t^5+56*n^2*p^3*t^3+24*n^2*p^2*t^4+106*n*p^4*t^3+192*n*p^3*t^4+90*
p^4*t^4-60*p^3*t^5-24*n^2*p^2*t^3-212*n*p^3*t^3-96*n*p^2*t^4-90*p^4*t^3-180*p^3*t^4-12*n^2*p^2*t^2-4*n^2*p*t^3+88*n*p^2*t^3+180*p^3*t^3+90*p^2*t^4+6*n^2*p^2*t+12*n^2*p*t^2+46*n*p^2*t^2+18*n*p*t^3-70*p^2*t^3-6*n^2*p*t-2*n^2*t^2-20*n*p^2*t-46*n*p*t^2-40*p^2*t^2-20*p*t^3+3*n^
2*t+20*n*p*t+7*n*t^2+14*p^2*t+40*p*t^2-n^2-10*n*t-14*p*t-5*t^2+3*n+7*t-2)/(n-1)/(4*n*p^2*t^2-4*n*p*t^2-10*p^2*t^2+10*p*t^2+2*n*t-n-5*t+2)/(p*t-1)/(p*t-t+1), p*t^2*(p-1)*(32*n^2*p^4*t^4-64*n^2*p^3*t^4-128*n*p^4*t^4+48*n^2*p^2*t^4+256*n*p^3*t^4+120*p^4*t^4-16*n^2*p*t^4-192*n
*p^2*t^4-240*p^3*t^4-12*n^2*p^2*t^2+64*n*p*t^4+180*p^2*t^4+12*n^2*p*t^2+8*n^2*t^3+42*n*p^2*t^2-60*p*t^4-12*n^2*t^2-42*n*p*t^2-32*n*t^3-30*p^2*t^2+2*n^2*t+46*n*t^2+30*p*t^2+30*t^3+n^2-7*n*t-40*t^2-3*n+5*t+2)/(n-1)/(4*n*p^2*t^2-4*n*p*t^2-10*p^2*t^2+10*p*t^2+2*n*t-n-5*t+2)/(p
*t-1)/(p*t-t+1), -2*p^2*t^4*(p-1)^2*(2*n-5)*(4*n*p^2*t^2-4*n*p*t^2-6*p^2*t^2+6*p*t^2+2*n*t-n-3*t+1)/(n-1)/(4*n*p^2*t^2-4*n*p*t^2-10*p^2*t^2+10*p*t^2+2*n*t-n-5*t+2)/(p*t-1)/(p*t-t+1)]:

if not (type(n1,integer) and n1>0) then
RETURN(FAIL):
fi:

if n1>=1 and n1<=3 then
L1[n1]:
else
normal(add(FtDfast(p,t,n1-i1)*subs(n=n1,A[i1]),i1=1..3)):
fi:

end:




#TFPstoryE(p,K): inputs a numeric or symbolic prob. p and integer K1 outputs an article about the quantity
#Expected family size in Talmudic Family Planning until reaching B boys AND G girls, it gives the
#two pure recurrences and a K by K matrix (as a list of lists) such that M[B][G] is the expected
#family size if you goal is to have B boys and G girls. Try:
#TFPstoryE(p,10);
TFPstoryE:=proc(p,K) local h,t,ope,H,T,i1,j1,L,t0:
t0:=time():
ope:=OPEAveF(h,t,H,T,p):

print(`A fast way to compute the Expected Number of Coin Tosses until you reach H heads AND T tails`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: You toss a coin whose probabilty of Heads is p and probability of Tails is 1-p and stop as soon as you have reached your goal of`):
print(`having at least h Heads AND t tails. Let L(h,t) be the expected number of tosses. Then L(h,t) satisfy the following pure recurrences`):
print(``):
print(``):print(L(h,t)=add(coeff(ope[1],T,-i1)*L(h,t-i1),i1=1..3)):
print(``):
print(``):print(L(h,t)=add(coeff(ope[2],H,-i1)* L(h-i1,t),i1=1..3)):
print(``):
print(`subject to the initial conditions `):
print(``):
print(seq(seq(L(i1,j1)=FaveF(i1,j1,p),j1=1..3),i1=1..3)):
print(``):
print(`and in Maple notation`):
print(``):
print(``):lprint(L(h,t)=add(coeff(ope[1],T,-i1)*L(h,t-i1),i1=1..3)):
print(``):
print(``):lprint(L(h,t)=add(coeff(ope[2],H,-i1)* L(h-i1,t),i1=1..3)):
print(``):
print(``):
print(seq(seq(L(i1,j1)=FaveF(i1,j1,p),j1=1..3),i1=1..3)):
print(``):

print(`These recurrences enable very fast compuations of the expected number of coin tosses until reahing the desired goal`):
print(``):
print(`For a fair coin, L(h,h)/h for multiples of 50 from 1 to `, 50*K, ` are `):
print(``):
print(seq(L(50*i1,50*i1)/(50*i1)=evalf(FaveF(50*i1,50*i1,1/2)/(50*i1)),i1=1..K)):

print(`If the prob. of a Head is 1/3 then, L(20*i,30*i)/(20*i) for i from 1 to`, K, ` are: `):
print(``):
print(seq(L(20*i1,40*i1)/(20*i1)=evalf(FaveF(20*i1,40*i1,1/3)/(20*i1)),i1=1..K)):
print(``):
print(`If the prob. of a Head is 2/5 then, L(20*i,30*i)/(20*i) for i from 1 to`, K, ` are: `):
print(``):
print(seq(L(20*i1,30*i1)/(20*i1)=evalf(FaveF(20*i1,30*i1,2/5)/(20*i1)),i1=1..K)):
print(``):
print(`-------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `seconds to produce `):
print(``):

end:



#TFPstoryGFd(p,t,K1,K2): inputs a numeric or symbolic prob. p a symbol t, and integer K outputs an article about the 
#prob. generating functions for the famiy sizes until reaching n boys AND n girls, it gives the
#recurrence (the same) for the case where you quit as soon as you have n babies of one gender and
#the case where you need at least n babies from both genders.
#then it compares the statistics of these two kinds of goals
#TFPstoryGFd(p,t,10,4);
TFPstoryGFd:=proc(p,t,K1,K2) local n,ope,N,i1,lu1,lu2,t0:
t0:=time():
ope:=OpeFt(p,t,n,N):
print(`A fast way to compute the Probability Generating Function for the Number of Coin Tosses until you reach n heads OR n tails and`):
print(`A fast way to compute the Probability Generating Function for the Number of Coin Tosses until you reach n heads AND n tails and`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: You toss a coin whose probabilty of Heads is p and probability of Tails is 1-p `):
print(`Let L1(n,t) be the probability generating function of the random variable: Number of coin tosses until you either get n heads or n tails`):
print(``):
print(`Let L2(n,t) be the probability generating function of the random variable: Number of coin tosses until you either get n heads AND n tails`):
print(``):
print(`Then surpisingly both these sequences of functions (the first polynomial, the second rational cunction, satisfies the SAME third-order recurrence`):
print(``):
print(``):print(add(coeff(ope,N,i1)*L1(n,t+i1),i1=0..3)=0):
print(``):
print(``):print(add(coeff(ope,N,i1)*L2(n,t+i1),i1=0..3)=0):
print(``):
print(`but with different initial conditions, of course.`):
print(``):
print(seq(L1(i1,t)=ftDfast(p,t,i1),i1=1..3)):
print(``):
print(seq(L2(i1,t)=FtDfast(p,t,i1),i1=1..3)):
print(``):
print(`In Maple notation: `):
print(``):
lprint(add(coeff(ope,N,i1)*L1(n,t+i1),i1=0..3)=0):
print(``):
lprint(add(coeff(ope,N,i1)*L2(n,t+i1),i1=0..3)=0):
print(``):
lprint(seq(L1(i1,t)=ftDfast(p,t,i1),i1=1..3)):
print(``):
lprint(seq(L2(i1,t)=FtDfast(p,t,i1),i1=1..3)):
print(``):


print(`These recurrences enable very fast compuations of these probabiltiy generating function, that enable us to compute expectation, variance and higher scaled moments. Let's assume a fair cooin.`):

for i1 from 1 to K1 do
print(``):
print(`-----------------------------------------------------`):
print(``):
print(`If your goal is`, 100*i1, `heads OR `, 100*i1, `tails then the expected duration, variance, up to the`, K2, `scaled-moments are `):
lu1:=ftDfast(1/2,t,i1*100):
print(``):
print(evalf(Stat(lu1,t,K2))):
print(``):
print(`If your goal is`, 100*i1, `heads AND `, 100*i1, `tails then the expected duration, variance, up to the`, K2, `scaled-moments are `):
lu2:=FtDfast(1/2,t,i1*100):
print(``):
print(evalf(Stat(lu2,t,K2))):
print(``):
od:
print(``):
print(`-------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `seconds to produce `):
print(``):


end:

#OPEx1x2(p,x1,x2,n,m,N,M): The pair of pure operators in n and m (with shift operators N and M) annihilating both
#f(n,m,p,1-p,x1,x2) and F(n,m,p,1-p,x1,x2) that enable fast computations. Try:
#OPEx1x2(p,x1,x2,n,m,N,M);
OPEx1x2:=proc(p,x1,x2,n,m,N,M)
[-p^2*x1^2*(n+m+1)/(p*x2-x2+1)/(n+2)+p*x1*(n*p*x1+n*p*x2-n*x2+(m+1)*p*x1+(m+1)*p*x2+2*n
-(m+1)*x2+m+3)/(p*x2-x2+1)/(n+2)*N-(n*p^2*x1*x2-n*p*x1*x2+(m+1)*p^2*x1*x2+2*n*p*x1+n*p*x2-(m+1)*p*x1*x2-n*x2+(m+3)*p*x1+2*p*x2+n-2*x2+2)/(p*x2-x2+1)/(n+2)*N^2+N^3,
x2^2*(p-1)^2*(m+n+1)/(p*x1-1)/(m+2)-x2*(p-1)*(m*p*x1+m*p*x2-m*x2+(n+1)*p*x1+(n+1)*p*x2-2*m-(n+1)*x2-n-3)/(p*x1-1)/(m+2)*M
+(m*p^2*x1*x2-m*p*x1*x2+(n+1)*p^2*x1*x2-m*p*x1-2*m*p*x2-(n+1)*p*x1*x2+2*m*x2-2*p*x1-(n+3)*p*x2+m+(n+3)*x2+2)/(p*x1-1)/(m+2)*M^2+M^3]:
end:






#Ffast(n1,m1,p,x1,x2): A fast way, using the third-order recurrences for  F(n1,m1,p,1-p,x1,x2). Try:
#F(6,3,p,1-p,x1,x2);
Ffast:=proc(n1,m1,p,x1,x2) local n,m,i1,L,kaM,kaN:
option remember:

kaN:=
[(m*p^2*x1*x2+n*p^2*x1*x2-m*p*x1*x2-n*p*x1*x2-2*p^2*x1*x2+m*p*x1+2*n*p*x1+n*p*x2+2*p*x1*x2-n*x2-3*p*x1-p*x2+n+x2-1)/(p*x2-x2+1)/(n-1),
-p*x1*(m*p*x1+m*p*x2+n*p*x1+n*p*x2-m*x2-n*x2-2*p*x1-2*p*x2+m+2*n+2*x2-3)/(p*x2-x2+1)/(n-1),
p^2*x1^2*(n-2+m)/(p*x2-x2+1)/(n-1)]:



kaM:=
[-(m*p^2*x1*x2+n*p^2*x1*x2-m*p*x1*x2-n*p*x1*x2-2*p^2*x1*x2-m*p*x1-2*m*p*x2-n*p*x2+2*p*x1*x2+2*m*x2+n*x2+p*x1+3*p*x2+m-3*x2-1)/(p*x1-1)/(m-1),
x2*(p-1)*(m*p*x1+m*p*x2+n*p*x1+n*p*x2-m*x2-n*x2-2*p*x1-2*p*x2-2*m-n+2*x2+3)/(p*x1-1)/(m-1),
-x2^2*(p-1)^2*(n-2+m)/(p*x1-1)/(m-1)]:

L:=
[[-(-1+p)*x2*p*x1*(p*x1-p*x2+x2-2)/(p*x2-x2+1)/(p*x1-1), (-1+p)^2*x2^2*p*x1*(p^2*x1^2-p^2*x1*x2+p*x1*x2-3*p*x1+2*p*x2-2*x2+3)/(p*x2-x2+1)/(p*x1-1)^2, -(-1+p)^3*x2^3*p*x1*(p^3*x1^3-p^3*x1^2*x2+p^2*x1^2*x2-4*p^2*x1^2+3*p^2*x1*x2-3*p*x1*x2+6*p*x1-3*p*x2+3*x2-4)/(p*x2-x2+1)/(p
*x1-1)^3], [-p^2*x1^2*(-1+p)*x2*(p^2*x1*x2-p^2*x2^2-p*x1*x2+2*p*x2^2+2*p*x1-3*p*x2-x2^2+3*x2-3)/(p*x2-x2+1)^2/(p*x1-1), (-1+p)^2*x2^2*p^2*x1^2*(2*p^3*x1^2*x2-2*p^3*x1*x2^2-2*p^2*x1^2*x2+4*p^2*x1*x2^2+3*p^2*x1^2-8*p^2*x1*x2+3*p^2*x2^2-2*p*x1*x2^2+8*p*x1*x2-6*p*x2^2-8*p*x1+8
*p*x2+3*x2^2-8*x2+6)/(p*x2-x2+1)^2/(p*x1-1)^2, -(-1+p)^3*x2^3*p^2*x1^2*(3*p^4*x1^3*x2-3*p^4*x1^2*x2^2-3*p^3*x1^3*x2+6*p^3*x1^2*x2^2+4*p^3*x1^3-15*p^3*x1^2*x2+8*p^3*x1*x2^2-3*p^2*x1^2*x2^2+15*p^2*x1^2*x2-16*p^2*x1*x2^2-15*p^2*x1^2+25*p^2*x1*x2-6*p^2*x2^2+8*p*x1*x2^2-25*p*x1
*x2+12*p*x2^2+20*p*x1-15*p*x2-6*x2^2+15*x2-10)/(p*x2-x2+1)^2/(p*x1-1)^3], [-p^3*x1^3*(-1+p)*x2*(p^3*x1*x2^2-p^3*x2^3-2*p^2*x1*x2^2+3*p^2*x2^3+3*p^2*x1*x2-4*p^2*x2^2+p*x1*x2^2-3*p*x2^3-3*p*x1*x2+8*p*x2^2+x2^3+3*p*x1-6*p*x2-4*x2^2+6*x2-4)/(p*x2-x2+1)^3/(p*x1-1), (-1+p)^2*x2^
2*p^3*x1^3*(3*p^4*x1^2*x2^2-3*p^4*x1*x2^3-6*p^3*x1^2*x2^2+9*p^3*x1*x2^3+8*p^3*x1^2*x2-15*p^3*x1*x2^2+4*p^3*x2^3+3*p^2*x1^2*x2^2-9*p^2*x1*x2^3-8*p^2*x1^2*x2+30*p^2*x1*x2^2-12*p^2*x2^3+3*p*x1*x2^3+6*p^2*x1^2-25*p^2*x1*x2+15*p^2*x2^2-15*p*x1*x2^2+12*p*x2^3+25*p*x1*x2-30*p*x2^
2-4*x2^3-15*p*x1+20*p*x2+15*x2^2-20*x2+10)/(p*x2-x2+1)^3/(p*x1-1)^2, -(-1+p)^3*x2^3*p^3*x1^3*(6*p^5*x1^3*x2^2-6*p^5*x1^2*x2^3-12*p^4*x1^3*x2^2+18*p^4*x1^2*x2^3+15*p^4*x1^3*x2-36*p^4*x1^2*x2^2+15*p^4*x1*x2^3+6*p^3*x1^3*x2^2-18*p^3*x1^2*x2^3-15*p^3*x1^3*x2+72*p^3*x1^2*x2^2-\
45*p^3*x1*x2^3+6*p^2*x1^2*x2^3+10*p^3*x1^3-63*p^3*x1^2*x2+63*p^3*x1*x2^2-10*p^3*x2^3-36*p^2*x1^2*x2^2+45*p^2*x1*x2^3+63*p^2*x1^2*x2-126*p^2*x1*x2^2+30*p^2*x2^3-15*p*x1*x2^3-36*p^2*x1^2+90*p^2*x1*x2-36*p^2*x2^2+63*p*x1*x2^2-30*p*x2^3-90*p*x1*x2+72*p*x2^2+10*x2^3+45*p*x1-45*
p*x2-36*x2^2+45*x2-20)/(p*x2-x2+1)^3/(p*x1-1)^3]]:

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

if n1<=3 and m1>=3  then
RETURN(normal(add(subs({n=n1,m=m1},kaM[i1])*Ffast(n1,m1-i1,p,x1,x2),i1=1..3))):
fi:


if n1>3  then
RETURN(normal(add(subs({n=n1,m=m1},kaN[i1])*Ffast(n1-i1,m1,p,x1,x2),i1=1..3))):
fi:

FAIL:
end:



#ffast(n1,m1,p,x1,x2): A fast way, using the third-order recurrences for  f(n1,m1,p,1-p,x1,x2). Try:
#f(6,3,p,1-p,x1,x2);
ffast:=proc(n1,m1,p,x1,x2) local n,m,i1,L,kaM,kaN:
option remember:

kaN:=
[(m*p^2*x1*x2+n*p^2*x1*x2-m*p*x1*x2-n*p*x1*x2-2*p^2*x1*x2+m*p*x1+2*n*p*x1+n*p*x2+2*p*x1*x2-n*x2-3*p*x1-p*x2+n+x2-1)/(p*x2-x2+1)/(n-1),
-p*x1*(m*p*x1+m*p*x2+n*p*x1+n*p*x2-m*x2-n*x2-2*p*x1-2*p*x2+m+2*n+2*x2-3)/(p*x2-x2+1)/(n-1),
p^2*x1^2*(n-2+m)/(p*x2-x2+1)/(n-1)]:



kaM:=
[-(m*p^2*x1*x2+n*p^2*x1*x2-m*p*x1*x2-n*p*x1*x2-2*p^2*x1*x2-m*p*x1-2*m*p*x2-n*p*x2+2*p*x1*x2+2*m*x2+n*x2+p*x1+3*p*x2+m-3*x2-1)/(p*x1-1)/(m-1),
x2*(p-1)*(m*p*x1+m*p*x2+n*p*x1+n*p*x2-m*x2-n*x2-2*p*x1-2*p*x2-2*m-n+2*x2+3)/(p*x1-1)/(m-1),
-x2^2*(p-1)^2*(n-2+m)/(p*x1-1)/(m-1)]:

L:=
[[(1-p)*x2+p*x1, (1-p)^2*x2^2+p*x1+p*x1*(1-p)*x2, (1-p)^3*x2^3+p*x1+p*x1*(1-p)*x2+p*x1*(1-p)^2*x2^2], [(1-p)*x2+p*x1*(1-p)*x2+p^2*x1^2, (1-p)^2*x2^2+2*p*x1*(1-p)^2*x2^2+p^2*x1^2+2*p^2*x1^2*(1-p)*x2, (1-p)^3*x2^3+3*p*x1*(1-p)^3*x2^3+p^2*x1^2+2*p^2*x1^2*(1-p)*x2+3*p^2*x1^2*(
1-p)^2*x2^2], [(1-p)*x2+p*x1*(1-p)*x2+p^2*x1^2*(1-p)*x2+p^3*x1^3, (1-p)^2*x2^2+2*p*x1*(1-p)^2*x2^2+3*p^2*x1^2*(1-p)^2*x2^2+p^3*x1^3+3*p^3*x1^3*(1-p)*x2, (1-p)^3*x2^3+3*p*x1*(1-p)^3*x2^3+6*p^2*x1^2*(1-p)^3*x2^3+p^3*x1^3+3*p^3*x1^3*(1-p)*x2+6*p^3*x1^3*(1-p)^2*x2^2]]:

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

if n1<=3 and m1>=3  then
RETURN(normal(add(subs({n=n1,m=m1},kaM[i1])*ffast(n1,m1-i1,p,x1,x2),i1=1..3))):
fi:


if n1>3  then
RETURN(normal(add(subs({n=n1,m=m1},kaN[i1])*ffast(n1-i1,m1,p,x1,x2),i1=1..3))):
fi:

FAIL:
end:



#faveF(n1,m1,p): A fast way to do fave(n1,m1,p), using the recurrence
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:=[[1, -p+2, p^2-3*p+3], [p+1, -2*p^2+2*p+2, 3*p^3-7*p^2+3*p+3], [p^2+p+1, -3*p^3+2*p^2+2*p+2, 6*p^4-12*p^3+3*p^2+3*p+3]]:

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

if n1<=3 and m1>=3  then
RETURN(normal(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(normal(add(subs({n=n1,m=m1},coeff(ope[2],N,-i1))*faveF(n1-i1,m1,p),i1=1..3))):
fi:

FAIL:
end:


#AveAndMoms2(f,x1,x2,N): Given a specific
#probability generating function
#for a bivariate distribution
#f(x1,x2) (a polynomial, rational function and possibly beyond)
#returns the averages of the first and second r.v. followed
#by a list-of-lists L
#such that L[r+1][s+1] is the (r,s)-mixed moment about the mean
#For example, try:
#AveAndMoms2(((1+x1+x2)/3)^100,x1,x2,4);

AveAndMoms2:=proc(f,x1,x2,N) local mu,F,memu1,memu2,i1,i2,T,F1,F11:
mu:=simplify(subs({x1=1,x2=1},f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs({x1=1,x2=1},diff(F,x1))):
memu2:=simplify(subs({x1=1,x2=1},diff(F,x2))):


F:=F/x1^memu1/x2^memu2:

F1:=F:
for i1 from 0 to N do
F11:=F1:
 for i2 from 0 to N do
  T[i1,i2]:=subs({x1=1,x2=1},F11):
  F11:=normal(x2*diff(F11,x2)):
 od:
 F1:=normal(x1*diff(F1,x1)):
od:

[memu1,memu2,[seq([seq(T[i1,i2],i2=0..N)],i1=0..N)]]:

end:


#AveAndMoms2S(f,x1,x2,N): Given a specific
#probability generating function
#for a bivariate distribution
#f(x1,x2) (a polynomial, rational function and possibly beyond)
#returns the averages of the first and second r.v. followed
#by a list-of-lists L
#such that L[r+1][s+1] is the SCALED (r,s)-mixed moment about the mean
#For example, try:
#AveAndMoms2S(((1+x1+x2)/3)^100,x1,x2,4);

AveAndMoms2S:=proc(f,x1,x2,N) local gu,gu1,gu2,gu3,i1,j1,v1,v2:

gu:=AveAndMoms2(f,x1,x2,N) :

gu1:=gu[1]:
gu2:=gu[2]:
gu3:=gu[3]:

v1:=gu3[1][3]:
v2:=gu3[3][1]:

gu3:=[seq([seq(gu3[i1+1][j1+1]/v1^(i1/2)/v2^(j1/2),j1=0..nops(gu3[i1+1])-1)],i1=0..nops(gu3)-1)]:

[[gu1,v1],[gu2,v2],gu3]:
end:




#TFPstoryGFbi(p,K)): inputs a numeric or symbolic prob. p and integer K outputs an article about the bi-variate
#prob. generating function for the number of heads and the number of tails until reaching either
#h Heads OR t TAILS
#and also when the goal is
#h Heads AND t TAILS
#where Pr(H)=p and Pr(T)=1-p
#It gives the linear recurrences (the SAME) and the respective initial conditions and
#then does some statistical analysis. Try:
#TFPstoryGFbi(p,100);
TFPstoryGFbi:=proc(p,K) local h,t,ope, i,i1,j1,L,t0,x1,x2,gu,lu,M,N:
t0:=time():

ope:=OPEx1x2(p,x1,x2,h,t,N,M):


print(`A fast way to compute the Bi-Variate Probability Generating Function for the Number of Heads and Tails  Coin Tosses until you reach H heads AND T tails (and also the case where you reach H heads OR T tails)`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: You toss a coin whose probabilty of Heads is p and probability of Tails is 1-p and stop as soon as you have reached your goal of`):
print(`having at least h Heads AND t tails. Let L1(h,t)=L1(h,t,p,x1,x2) be the generating function whoere coefficient of x1^h1*x2^t1 is the exact probability`):
print(`that when you have reached your goal of at least h Heads AND t Tails (of course either h1=h or t1=t)`):
print(`Analogously, let L2(h,t)=L2(h,t,p,x1,x2) be the bi-variate  prob. generating function where you stop as soon as you have reached`):
print(`h Heads OR t Tails `):
print(``):
print(`Then BOTH of these discrete functions of h and t satisfy the SAME third-order recurrences`):

print(``):
print(add(coeff(ope[1],N,i1)*L(h+i1,t),i1=0..degree(ope[1],N))=0):
print(``):
print(add(coeff(ope[2],M,i1)*L(h,t+i1),i1=0..degree(ope[2],M))=0):
print(``):
print(`subject to the initial conditions, respectively `):
print(``):
print(seq(seq(L1(i1,j1)=Ffast(i1,j1,p,x1,x2),j1=1..3),i1=1..3)):
print(``):
print(``):
print(seq(seq(L2(i1,j1)=ffast(i1,j1,p,x1,x2),j1=1..3),i1=1..3)):
print(``):

print(`and in Maple notation`):


print(``):
lprint(add(coeff(ope[1],N,i1)*L(h+i1,t),i1=0..degree(ope[1],N))=0):
print(``):
lprint(add(coeff(ope[2],M,i1)*L(h,t+i1),i1=0..degree(ope[2],M))=0):
print(``):
print(`subject to the initial conditions, respectively `):
print(``):
lprint(seq(seq(L1(i1,j1)=Ffast(i1,j1,p,x1,x2),j1=1..3),i1=1..3)):
print(``):
print(``):
lprint(seq(seq(L2(i1,j1)=ffast(i1,j1,p,x1,x2),j1=1..3),i1=1..3)):
print(``):



print(`These recurrences enable very fast compuations of the expected number of these bi-variate prob. gen. functions and deducing relevant statistics.`):

print(` Chapter I`):
print(``):
print(`Let's illustrate it with p=1/2 and (m,n)=(50*i,50*i) for i from 1 to`, K):

for i from 1 to K do

print(`if the goal is `, 50*i , `Heads OR `, 50*i , `Tails (you are happy once you have one of them) `):

gu:=ffast(50*i,50*i,1/2,x1,x2):

lu:=evalf(AveAndMoms2S(gu,x1,x2,4)):

print(`The expected number of Heads is`,lu[1][1], `and the standard-deviation is `, lu[1][2]):
print(``):
print(`The expected number of Tails is`,lu[2][1], `and the standard-deviation is `, lu[2][2]):
print(``):
print(`of course they are the same by symmetry `):
print(``):
print(`The skewness and kurtosis of the number of Heads are`, lu[3][1][4], lu[3][1][5], ` respectively `):
print(``):
print(`The skewness and kurtosis of the number of Tails are`, lu[3][4][1], lu[3][5][1], `respectively `):
print(``):
print(`of course they are the same by symmetry `):
print(``):
print(`More interesting is the coefficient of correlation that is `, lu[3][2][2]):
print(``):

print(`On the other hand if the goal is (at least) `, 50*i , `Heads AND  (at least)`, 50*i , `Tails (you are only happy once you have BOTH  goals `):

gu:=Ffast(50*i,50*i,1/2,x1,x2):

lu:=evalf(AveAndMoms2S(gu,x1,x2,4)):

print(`The expected number of Heads is`,lu[1][1], `and the standard-deviation is `, lu[1][2]):
print(``):
print(`The expected number of Tails is`,lu[2][1], `and the standard-deviation is `, lu[2][2]):
print(``):
print(`of course they are the same by symmetry `):
print(``):
print(`The skewness and kurtosis of the number of Heads are`, lu[3][1][4], lu[3][1][5], ` respectively `):
print(``):
print(`The skewness and kurtosis of the number of Tails are`, lu[3][4][1], lu[3][5][1], `respectively `):
print(``):
print(`of course they are the same by symmetry `):
print(``):
print(`More interesting is the coefficient of correlation that is `, lu[3][2][2]):
print(``):
od:

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


print(` Chapter II`):
print(``):
print(`Let's illustrate it with p=1/3 and (m,n)=(20*i,40*i) for i from 1 to`, K):

for i from 1 to K do

print(`if the goal is `, 20*i , `Heads OR `, 40*i , `Tails (you are happy once you have one of them) `):

gu:=ffast(20*i,40*i,1/3,x1,x2):

lu:=evalf(AveAndMoms2S(gu,x1,x2,4)):

print(`The expected number of Heads is`,lu[1][1], `and the standard-deviation is `, lu[1][2]):
print(``):
print(`The expected number of Tails is`,lu[2][1], `and the standard-deviation is `, lu[2][2]):
print(``):
print(`Note that the ratio of the expectations is`, lu[2][1]/lu[1][1]):
print(``):
print(`The skewness and kurtosis of the number of Heads are`, lu[3][1][4], lu[3][1][5], ` respectively `):
print(``):
print(`The skewness and kurtosis of the number of Tails are`, lu[3][4][1], lu[3][5][1], `respectively `):
print(``):
print(``):
print(`More interesting is the coefficient of correlation that is `, lu[3][2][2]):
print(``):

print(`On the other hand if the goal is (at least) `, 20*i , `Heads AND  (at least)`, 40*i , `Tails (you are only happy once you have BOTH  goals `):

gu:=Ffast(20*i,40*i,1/3,x1,x2):

lu:=evalf(AveAndMoms2S(gu,x1,x2,4)):

print(`The expected number of Heads is`,lu[1][1], `and the  standard-deviation is `, lu[1][2]):
print(``):
print(`The expected number of Tails is`,lu[2][1], `and the standard-deviation is `, lu[2][2]):
print(``):
print(`Note that the ratio of the expectations is`, lu[2][1]/lu[1][1]):
print(``):
print(`The skewness and kurtosis of the number of Heads are`, lu[3][1][4], lu[3][1][5], ` respectively `):
print(``):
print(`The skewness and kurtosis of the number of Tails are`, lu[3][4][1], lu[3][5][1], `respectively `):
print(``):
print(``):
print(`More interesting is the coefficient of correlation that is `, lu[3][2][2]):
print(``):
od:

print(``):


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


print(` Chapter III`):
print(``):
print(`Let's illustrate it with p=1/4 and (m,n)=(20*i,60*i) for i from 1 to`, K):

for i from 1 to K do

print(`if the goal is `, 20*i , `Heads OR `, 60*i , `Tails (you are happy once you have one of them) `):

gu:=ffast(20*i,60*i,1/4,x1,x2):

lu:=evalf(AveAndMoms2S(gu,x1,x2,4)):

print(`The expected number of Heads is`,lu[1][1], `and the standard-deviation is `, lu[1][2]):
print(``):
print(`The expected number of Tails is`,lu[2][1], `and the standard-deviation is `, lu[2][2]):
print(``):
print(`Note that the ratio of the expectations is`, lu[2][1]/lu[1][1]):
print(``):
print(`The skewness and kurtosis of the number of Heads are`, lu[3][1][4], lu[3][1][5], ` respectively `):
print(``):
print(`The skewness and kurtosis of the number of Tails are`, lu[3][4][1], lu[3][5][1], `respectively `):
print(``):
print(``):
print(`More interesting is the coefficient of correlation that is `, lu[3][2][2]):
print(``):

print(`On the other hand if the goal is (at least) `, 20*i , `Heads AND  (at least)`, 40*i , `Tails (you are only happy once you have BOTH  goals `):

gu:=Ffast(20*i,60*i,1/4,x1,x2):

lu:=evalf(AveAndMoms2S(gu,x1,x2,4)):

print(`The expected number of Heads is`,lu[1][1], `and the  standard-deviation is `, lu[1][2]):
print(``):
print(`The expected number of Tails is`,lu[2][1], `and the standard-deviation is `, lu[2][2]):
print(``):
print(`Note that the ratio of the expectations is`, lu[2][1]/lu[1][1]):
print(``):
print(`The skewness and kurtosis of the number of Heads are`, lu[3][1][4], lu[3][1][5], ` respectively `):
print(``):
print(`The skewness and kurtosis of the number of Tails are`, lu[3][4][1], lu[3][5][1], `respectively `):
print(``):
print(``):
print(`More interesting is the coefficient of correlation that is `, lu[3][2][2]):
print(``):
od:

print(``):


print(` Chapter IV`):
print(``):
print(`Let's illustrate it with p=2/5 and (m,n)=(20*i,30*i) for i from 1 to`, K):

for i from 1 to K do

print(`if the goal is `, 20*i , `Heads OR `, 30*i , `Tails (you are happy once you have one of them) `):

gu:=ffast(20*i,30*i,2/5,x1,x2):

lu:=evalf(AveAndMoms2S(gu,x1,x2,4)):

print(`The expected number of Heads is`,lu[1][1], `and the standard-deviation is `, lu[1][2]):
print(``):
print(`The expected number of Tails is`,lu[2][1], `and the standard-deviation is `, lu[2][2]):
print(``):
print(`Note that the ratio of the expectations is`, lu[2][1]/lu[1][1]):
print(``):
print(`The skewness and kurtosis of the number of Heads are`, lu[3][1][4], lu[3][1][5], ` respectively `):
print(``):
print(`The skewness and kurtosis of the number of Tails are`, lu[3][4][1], lu[3][5][1], `respectively `):
print(``):
print(``):
print(`More interesting is the coefficient of correlation that is `, lu[3][2][2]):
print(``):

print(`On the other hand if the goal is (at least) `, 20*i , `Heads AND  (at least)`, 30*i , `Tails (you are only happy once you have BOTH  goals `):

gu:=Ffast(20*i,30*i,2/5,x1,x2):

lu:=evalf(AveAndMoms2S(gu,x1,x2,4)):

print(`The expected number of Heads is`,lu[1][1], `and the  standard-deviation is `, lu[1][2]):
print(``):
print(`The expected number of Tails is`,lu[2][1], `and the standard-deviation is `, lu[2][2]):
print(``):
print(`Note that the ratio of the expectations is`, lu[2][1]/lu[1][1]):
print(``):
print(`The skewness and kurtosis of the number of Heads are`, lu[3][1][4], lu[3][1][5], ` respectively `):
print(``):
print(`The skewness and kurtosis of the number of Tails are`, lu[3][4][1], lu[3][5][1], `respectively `):
print(``):
print(``):
print(`More interesting is the coefficient of correlation that is `, lu[3][2][2]):
print(``):
od:

print(``):

print(`This ends this paper that took`, time()-t0, `seconds to generate. `):
print(``):

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

print(``):
print(`This took`, time()-t0, `seconds. `):

end:


#rf(a,n): The raising factorial a(a+1)...(a+n-1). Try:
#rf(1/2,7);
rf:=proc(a,n) local i: mul(a+i,i=0..n-1):end:

#CheckEB(a,b,R): Checks the explicit formulas for  Fave(a*n,b*n,a/(a+b)) and fave(a*n,b*n,a/(a+b)) for n from 1 to R. Try:
#CheckEB(3,11,20);
CheckEB:=proc(a,b,R) local n,L,L1,M,M1:

L:=[seq(FaveF(a*n,b*n,a/(a+b)),n=1..R)]:
#L1:=[seq((a+b)*n+(a+b)*mul(rf(i/(a+b),n),i=1..a+b-1)/(n-1)!/mul(rf(i/a,n),i=1..a-1)/mul(rf(i/b,n),i=1..b-1),n=1..R)]:
L1:=[seq( (a+b)*n+ (a+b)*n*binomial((a+b)*n,a*n)*((a^a*b^b)/(a+b)^(a+b))^n,n=1..R)]:
M:=[seq(faveF(a*n,b*n,a/(a+b)),n=1..R)]:
M1:=[seq( (a+b)*n- (a+b)*n*binomial((a+b)*n,a*n)*((a^a*b^b)/(a+b)^(a+b))^n,n=1..R)]:

evalb(L=L1) and evalb(M=M1):

end:








#FairMom(n,r): The r-th moment of the number of coin-tosses of a fair coin until reaching n Heads or n Tails. Try:
#FairMom(n,2);
FairMom:=proc(n,r) local b1:
2*sum((b1+n)^r*binomial(n+b1-1,n-1)/2^(b1+n),b1=0..n-1):
end:

#FairMomsC(n,R): The expectation, variance, and central moments for reaching n Heads for n Tails with a fair coin up to the R-th moment.
#It also returns the limits of the scaled moments and the floating point of the latter.Try:
#FairMomsC(n,4);
FairMomsC:=proc(n,R) local r,lu,mu,gu,ka,i1,vu:
lu:=[seq(FairMom(n,r),r=1..R)]:
mu:=lu[1]:
gu:=[mu]:
for r from 2 to R do
ka:=normal(add(lu[r-i1]*(-mu)^i1*binomial(r,i1),i1=0..r-1)+(-mu)^r):
gu:=[op(gu),ka]:
od:
gu:
vu:=[limit(sqrt(gu[2])/gu[1],n=infinity)]:
vu:=[op(vu),seq(limit(gu[i1]/gu[2]^(i1/2),n=infinity),i1=2..R)]:

[gu,simplify(vu),evalf(vu)]:

end:



#FairMomStoryOld(n,R): An article about the time it takes to reach n Heads or n Tails with a fair coin. Try:
#FairMomStoryOld(n,4);
FairMomStoryOld:=proc(n,R) local gu,t0,r:
gu:=FairMomsC(n,R):
t0:=time():
print(`The Expectation, Variance, and first`, R, `moments of the random variable: reaching n Heads or n Tails Tossing a Fair Coin`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(` You toss a fair coin and stop as soon as you have reached n Heads Or n Tails, How long should it take?`):
print(``):
print(`The expectation of this random variable is`):
print(``):
print(gu[1][1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1][1]):
print(``):
print(``):
print(`The variance of this random variable is`):
print(``):
print(``):
print(gu[1][2]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1][2]):
print(``):
print(`Hence the limit of the coefficient of variation is`, gu[2][1]):

for r from 3 to R do
print(`The `, r, `-th central limit is `):
print(``):
print(gu[1][r]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1][r]):
print(``):
print(`and the limit of the` , r, `scaled moment is `):
print(``):
print(gu[2][r]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[2][r]):
print(``):
print(`and in deminals`):
lprint(evalf(gu[2][r])):

od:
print(``):
print(`-----------------------------------`):
print(``):
print(`To sum up the first`,  R , ` moments are`):
print(``):
lprint(gu[1]):
print(``):
print(`The limits  of the scaled central moments, starting at the third  are:`):
print(``):
lprint([op(3..nops(gu[2]),gu[2])]):
print(``):
print(`and  in Floats it is`):
print(``):
lprint(evalf([op(3..nops(gu[2]),gu[2])])):
print(``):
print(`On the other hand `):
print(``):
print(`The scaled central moments of the negative absolute value of the Normal distribution whose pdf is exp(-x^2/2)/sqrt(Pi/2)), supported in [-infinity..0] are`):
print(``):
print(AbsNorMomsN(R)[2]):
print(``):
print(`Yea! These are the same!`):
print(``):
print(`Hence we have a purely elementary and discete confirmation that the scaled entral moments of our random variable "number of coin tosses" of a fair coin until reaching n Heads or n Tails  tends to those of minus the absolute value `):
print(`of the normal distribution, at least they match up to the first`, R, ` moments. `):
print(`-----------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `to produce `):

end:

#plotDist(f,x,K): Given a prob. gen. function f(x) that has a 
#Taylor series
#for a discrete r.v.
#plots its normalized version (X-mu)/sig between mu-K1*sig and
#mu+K2*sig
#For example, try:
#plotDist((1+x)^40,x,5,5);
plotDist:=proc(f,x,K1,K2) local mu,f1,lu,gu,sig,i,j1:
f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
gu:=f1/x^mu:
sig:=sqrt(subs(x=1,x*diff(x*diff(gu,x),x))):

lu:=taylor(f1,x=0,trunc(mu)+K2*trunc(sig)+10):

lu:=[seq([i,coeff(lu,x,i)],i=max(0,trunc(mu-K1*sig))..trunc(mu+K2*sig))]:

lu:=evalf([seq([(lu[j1][1]-mu)/sig,lu[j1][2]*sig],j1=1..nops(lu))]):

plot(lu):
#plot(lu,scaling=constrained):

end:

#AbsNorMoms(R): The mean, variance and the Scaled centered moments of |N(0,1)|.
#Try:
#AbsNorMoms(6);
AbsNorMoms:=proc(R) local gu,i,x,mu,lu,r,ka,i1,va:
lu:=[seq(int(x^i*exp(-(x^2/2))/sqrt(Pi/2),x=0..infinity),i=1..R)]:
mu:=lu[1]:
gu:=[mu]:
for r from 2 to R do
ka:=add(lu[r-i1]*(-mu)^i1*binomial(r,i1),i1=0..r-1)+(-mu)^r:
gu:=[op(gu),ka]:
od:
va:=gu[2]:
simplify(normal([[mu, va], [seq(gu[i]/va^(i/2),i=3..R)]])):

end:


#AbsNorMomsN(R): The mean, variance and the Scaled centered moments of -|N(0,1)|.
#Try:
#AbsNorMomsN(6);
AbsNorMomsN:=proc(R) local gu,i,x,mu,lu,r,ka,i1,va:
lu:=[seq(int(x^i*exp(-(x^2/2))/sqrt(Pi/2),x=-infinity..0),i=1..R)]:
mu:=lu[1]:
gu:=[mu]:
for r from 2 to R do
ka:=add(lu[r-i1]*(-mu)^i1*binomial(r,i1),i1=0..r-1)+(-mu)^r:
gu:=[op(gu),ka]:
od:
va:=gu[2]:
[[mu, va],[seq(simplify(gu[i]/va^(i/2)),i=3..R)]]:

end:




#FairFacMom(n,r): The r-th Factorial moment of the number of coin-tosses of a fair coin until reaching n Heads or n Tails. Try:
#FairFacMom(n,2);
FairFacMom:=proc(n,r) local b1:
2*sum(r!*binomial(b1+n,r)*binomial(n+b1-1,n-1)/2^(b1+n),b1=0..n-1):
end:


#FairFacMomF(n,r): The r-th Factorial moment of the number of coin-tosses of a fair coin until reaching n Heads or n Tails, done FAST (using the recurrence). Try:
#FairFacMomF(n,2);
FairFacMomF:=proc(n,r) local C:
option remember:

subs(C=binomial(2*n,n+1)*(n+1)/4^n, FairFacMomF1(n,r,C)):
#if r=0 then
#1:
#elif r=1 then
#2*n-2*(2*n)!/(n!*(n-1)!)/4^n:
#else
#simplify(2*n*FairFacMomF(n,r-1)+ (r-1)*(r-2)*FairFacMomF(n,r-2)-4*4^(-n)*n*(n+1)*binomial(2*n,n+1)*binomial(2*n-1,r-2)*(r-2)!):
#fi:

end:

#FairFacMomF1(n,r,C): The r-th Factorial moment of the number of coin-tosses of a fair coin until reaching n Heads or n Tails, done FAST (using the recurrence). Try:
#FairFacMomF1(n,2,C);
FairFacMomF1:=proc(n,r,C)
option remember:
#C=binomial(2*n,n+1)*(n+1)/4^n:
if r=0 then
1:
elif r=1 then
2*n-2*C
else
expand(2*n*FairFacMomF1(n,r-1,C)+ (r-1)*(r-2)*FairFacMomF1(n,r-2,C)-4*n*C*binomial(2*n-1,r-2)*(r-2)!):
fi:

end:



#FairMomF1(n,r,C): a fast way to compute FairMom(n,r) in terms of C:=binomial(2*n,n+1)*(n+1)/4^n. Try:
#FairMomF1(n,4,C);
FairMomF1:=proc(n,r,C) local i:
expand(add(stirling2(r,i)*FairFacMomF1(n,i,C),i=0..r)):
end:


#FairMomF(n,r): a fast way to compute FairMom(n,r) 
#FairMomF(n,4);
FairMomF:=proc(n,r) local C:

subs(C=binomial(2*n,n+1)*(n+1)/4^n,FairMomF1(n,r,C)):
end:




#FairMomsCF1(n,R,C): The expectation, variance, and central moments for reaching n Heads for n Tails with a fair coin up to the R-th moment.
#in terms of C:=binomial(2*n,n+1)*(n+1)/4^n:
#It also returns the limits of the scaled moments and the floating point of the latter.Try:
#FairMomsCF1(n,4,C);
FairMomsCF1:=proc(n,R,C) local r,lu,mu,gu,ka,i1,vu,vu1,i:
lu:=[seq(FairMomF1(n,r,C),r=1..R)]:
mu:=lu[1]:
gu:=[mu]:
for r from 2 to R do
ka:=normal(add(lu[r-i1]*(-mu)^i1*binomial(r,i1),i1=0..r-1)+(-mu)^r):
gu:=[op(gu),ka]:
od:
vu:=[sqrt(gu[2])/gu[1]]:


vu:=[op(vu),seq(normal(gu[i1]/gu[2]^(i1/2)),i1=2..R)]:

vu1:=normal(subs(C=sqrt(n)/sqrt(Pi),vu)):

vu1:=[seq(limit(vu1[i],n=infinity),i=3..nops(vu1))]:

[gu,simplify(vu1),evalf(vu1)]:

end:




#FairMomStory(n,R,S): An article about the time it takes to reach n Heads or n Tails with a fair coin.
#It gives the explicit expressions for the moments, and the scaled limits from 1 to R
#and for those from R+1 to S just checks that the limit of the scaled central moments are
#the same as those of -|N(0,1)| Try:
#FairMomStory(n,4,20);
FairMomStory:=proc(n,R,S) local gu,t0,r,C,C1:
C1:=binomial(2*n,n+1)*(n+1)/4^n:
gu:=FairMomsCF1(n,R,C):



t0:=time():
print(`The Expectation, Variance, and first`, R, `moments of the random variable: reaching n Heads or n Tails By Tossing a Fair Coin`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let C denote `, C1):
print(``):
print(`Note that it is asymptotically`, sqrt(n/Pi), `plus O(1/sqrt(n)) `):
print(``):

print(` You toss a fair coin and stop as soon as you have reached n Heads OR n Tails, How long should it take?`):
print(``):
print(`The expectation of this random variable is`):
print(``):
print(gu[1][1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1][1]):
print(``):
print(``):
print(`The variance of this random variable is`):
print(``):
print(``):
print(gu[1][2]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1][2]):
print(``):
print(`Hence the limit of the coefficient of variation is`, limit(sqrt(gu[1][2])/gu[1][1],n=infinity)):
print(``):
for r from 3 to R do
print(`The `, r, `-th central limit is , in terms of C= binomial(2*n,n+1)*(n+1)/4^n, is:`):
print(``):
print(gu[1][r]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1][r]):
print(``):
print(`and the limit of the` , r, `scaled moment is `):
print(``):
print(gu[2][r-2]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[2][r-2]):
print(``):
print(`and in deminals`):
lprint(evalf(gu[2][r-2])):

od:
print(``):
print(`-----------------------------------`):
print(``):
print(`To sum up the first`,  R , ` moments are`):
print(``):
lprint(gu[1]):
print(``):
print(`The limits  of the scaled central moments, starting at the third  are:`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu[2]):
print(``):
print(`and  in Floats it is`):
print(``):
lprint(evalf(gu[2])):
print(``):
print(`On the other hand `):
print(``):
print(`The scaled central moments of the negative absolute value of the Normal distribution whose pdf is exp(-x^2/2)/sqrt(Pi/2)), supported in [-infinity..0] are`):
print(``):

print(AbsNorMomsN(R)[2]):
print(``):
print(`Yea! These are the same!`):
print(``):
print(`It would be too complicated to actually list the explicit expressions for the moments from`, R+1 , ` to `, S, `but let's check that `):
print(`The limits  of the scaled central moments, up to the`, S):
print(`coincide with those  of the negative absolute value of the Normal distribution whose pdf is exp(-x^2/2)/sqrt(Pi/2)), supported in [-infinity..0] are`):
print(``):
print(`Let's see:`):
print(``):

if evalb(AbsNorMomsN(S)[2]=FairMomsCF1(n,S,C)[2]) then
print(``):
print(`Yea, they are the same!`):
print(``):

print(`Hence we have a purely elementary and discete confirmation that the scaled entral moments of our random variable "number of coin tosses" of a fair coin until reaching n Heads or n Tails  tends to those of minus the absolute value `):
print(`of the normal distribution, at least they match up to the first`, S, ` moments. `):

else
print(`Too bad, they are not`):
fi:

print(`-----------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `to produce `):

end: