######################################################################
## GDW.txt Save this file as  GDW.txt  to use it,                    #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read GWD.txt <Enter>                              #
## Then follow the instructions given there                          #
##                                                                   #
## Written by AJ Bu and Doron Zeilberger, Rutgers University ,       #
######################################################################

with(combinat):

print(`First Written: April, 2023: tested for Maple 20 `):
print():
print(`This is GWD.txt, A Maple package`):
print(`for the enumeration of Generalized Dyck Walks using Algebraic Schemes generated automatically`):

print(`accompanying AJ Bu and Doron Zeilberger's  article: `):
print(` Using Symbolic Computation to Explore  Generalized Dyck Paths and Their Areas`):

print():
print(`The most current version of this Maple package is available at:`):
print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/GWD.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`-------------------------------------------------------------------------------------------`):
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):

print(`-------------------------------------------------------------------------------------------`):
print(`For a list of the supporting functions type: ezra1();`):
print(`  For specific help type "ezra(procedure_name);" `):
print(`-------------------------------------------------------------------------------------------`):


print(`-------------------------------------------------------------------------------------------`):
print(`For a list of the  Area procedures type: ezraA();`):
print(`  For specific help type "ezra(procedure_name);" `):
print(`-------------------------------------------------------------------------------------------`):


print(`-------------------------------------------------------------------------------------------`):
print(`For a list of the  Support Area procedures type: ezraA1();`):
print(`  For specific help type "ezra(procedure_name);" `):
print(`-------------------------------------------------------------------------------------------`):



print(`-------------------------------------------------------------------------------------------`):
print(`For a list of the CHECKING functions type: ezraC();`):
print(`  For specific help type "ezra(procedure_name);" `):
print(`-------------------------------------------------------------------------------------------`):


print(`-------------------------------------------------------------------------------------------`):
print(`For a list of the STORY functions type: ezraSt();`):
print(`  For specific help type "ezra(procedure_name);" `):
print(`-------------------------------------------------------------------------------------------`):



print(`-------------------------------------------------------------------------------------------`):
print(`For a list of the SIMULATION functions type: ezraSi();`):
print(`  For specific help type "ezra(procedure_name);" `):
print(`-------------------------------------------------------------------------------------------`):



print():
 
#Digits:=100:




ezraA:=proc()
if args=NULL then

print(` The MAIN AREA procedures are: BookA, BookA2, qEqGFt, qEqGFtS, Paper1A, Paper11A, Paper2A  `):
print(`  `):
else
ezra(args):
fi:

end:


ezraA1:=proc()
if args=NULL then

print(` The SUPPORTING AREA procedures are: ModeP, qMakeEqFt,  qMakeEqGt, qMakeSyst, qSeqS `):
print(`  `):
else
ezra(args):
fi:

end:

ezraSi:=proc()
if args=NULL then

print(` The SIMULATION procedures are:`):
print(`  LD, OneTrip, SimuGF, SimuStat,  Stat  `):
else
ezra(args):
fi:

end:

ezraC:=proc()
if args=NULL then

print(` The Checking procedures are:`):
print(`  CheckEqGF,   CheckEqGFt, CheckEqP  `):
else
ezra(args):
fi:

end:

ezraSt:=proc()
if args=NULL then

print(` The STORY procedures are:`):
print(`  Book, Paper1, Paper11, Paper1p, Paper1pS, Paper2pS `):
else
ezra(args):
fi:

end:

ezra1:=proc()
if args=NULL then

print(` The supporting procedures are:`):
print(` AlgToSeq, AveSD, DecoW, DerEqP, DerEqPk, DerF, EqGFe, Findrec, GFpr, GFstx, IsStrict, MakeEqF,   MakeSys,  MakeSysE, MakeSysT, MakeSysT1, NuWabn, PWabn, PWE, SeqFromSysEq, WabnSdirect,  Wabn, WabnS, Wt `):
else
ezra(args):
fi:

end:

ezra:=proc()
if args=NULL then

print(` GWD.txt: A Maple package for  generating and processing Algebraic Schemes for Enumerating Generalized Dyck Paths`):

print(`The MAIN procedures are:   AveAndMomsP, AveAndMomsPc, EqGF, EqGFt, EqGFtS, EqP, SeqFromEq, SeqS, SeqSe, SeqSee `):

print(``):


elif nops([args])=1 and op(1,[args])=AlgToSeq then
print(`AlgToSeq(eq,X,t,K): Given an algebraic equation in terms of X=X(t) and t, finds the first K terms starting at n=0. Try:`):
print(`AlgToSeq(X-1-t*X^2,X,t,30);`):

elif nops([args])=1 and op(1,[args])=AveAndMomsP then
print(`AveAndMomsP(P,K): Inputs a losing (or symbolic) probability distribution P, finds a list of length K whose first entry is the expected duration until out of the game, the second is the standard deviation, followed by the 3rd through the`):
print(`K-th scaled moments. Try:`):
print(`AveAndMomsP([[1,2/5],[0,1/10],[-1,1/2]],4);`):


elif nops([args])=1 and op(1,[args])=AveAndMomsPc then
print(`AveAndMomsPc(st,P,K,L): Inputs a losing NUMERIC LOSING probability distribution P, a positive integer K, and a a (usually) large positive integer L`):
print(`STARTING with capital st`):
print(`outputs the a list of length K whose first entry is the CONDITIONAL expected duration until out of the game CONDITIONED that you exited after L steps`):
print(`the second is the CONDITIONAL standard deviation, followed by the 3rd through the`):
print(`K-th scaled moments. It also outputs the probability that you will exit in <=L rounds. Try:`):
print(`AveAndMomsPc(0,[[1,4/10],[0,1/10],[-1,1/2]],4,1000);`):


elif nops([args])=1 and op(1,[args])=AveSD then
print(`AveSD(P): The exact average and standard-deviation  of duration with NUMERIC prob. dist. P. Try:`):
print(`AveSD([[1,1/4],[0,1/4],[-1,1/2]]);`):
print(`AveSD([[2,1/8],[1,1/8],[0,1/8],[-1,5/16],[-2,5/16]]);`):

elif nops([args])=1 and op(1,[args])=Book then
print(`Book(M,MaxC): Inputs a positive integer M and outputs a book about enumearing generalized Dyck path with all possible alphabets consisting of both negative and positive integes and 0. It also`):
print(`outputs recurrences if their complexity is <=MaxC. Try:`):
print(`Book(2,10);`):

elif nops([args])=1 and op(1,[args])=BookA then
print(`BookA(M): Inputs a positive integer M and outputs a book about enumearing generalized Dyck path with all possible alphabets consisting of both negative and positive integes and 0.`):
print(`BookA(1);`):


elif nops([args])=1 and op(1,[args])=BookA2 then
print(`BookA2(M,MaxC): Inputs a positive integer M and outputs a book about enumearing generalized Dyck path as well as the sum of the areas and area-squared`):
print(`with all possible alphabets consisting of both negative and positive integes and 0. It also`):
print(`outputs recurrences if their complexity is <=MaxC. Try:`):
print(`BookA2(2,10);`):

elif nops([args])=1 and op(1,[args])=CheckEqGF then
print(`CheckEqGF(S,K): Checks procedure EqGF(S,X,x) up to K terms. Try:`):
print(`CheckEqGF({1,0,-1},10);`):


elif nops([args])=1 and op(1,[args])=CheckEqGFt then
print(`CheckEqGFt(S,K): Checks procedure EqGFt(S,X,t) up to K terms. Try:`):
print(`CheckEqGFt({1,0,-1},10);`):


elif nops([args])=1 and op(1,[args])=CheckEqGFt then
print(`CheckEqGFt(S,K): Checks procedure EqGFt(S,X,t) up to K terms. Try:`):
print(`CheckEqGFt({1,0,-1},10);`):

elif nops([args])=1 and op(1,[args])=CheckEqP then
print(`CheckEqP(P,K): Checks procedure EqP(P,X,t) up to K terms. Try:`):
print(`CheckEqP([[1,1/4],[0,1/4],[-1,1/2]],10);`):

elif nops([args])=1 and op(1,[args])=DecoW then
print(`DecoW(w,a): Given a walk that starts at [0,a] outputs [w1,w2] where w1 is the  largest prefix that is strictly above the x-axis, followed by the rest so w=w1w2. Try`):
print(` DecoW([1,-1,1,-1],0);`):


elif nops([args])=1 and op(1,[args])=DerEqP then
print(`DerEqP(P,X,t): an expression for t*X'(t) if X(t) is the prob. gen. function for the Genralized Dyck path process. Try:`):
print(`DerEqP([[1,a],[0,b],[-1,1-a-b]],X,t);`):

elif nops([args])=1 and op(1,[args])=DerEqPk then
print(`DerEqPk(P,X,t,k): an expression for (t*d/dt)^k X(t) if X(t) is the prob. gen. function for the Genralized Dyck path process. Try:`):
print(`DerEqPk([[1,a],[0,b],[-1,1-a-b]],X,t,2);`):

elif nops([args])=1 and op(1,[args])=DerF then
print(`DerF(F,X,t): inputs a polynomial f in X and t, uses implicit differention to finds an explicit expression for X'(t), in terms of X and t where X stands for X(t)`):
print(`Try: DerF(X-1-t*X^2,X,t);`):

elif nops([args])=1 and op(1,[args])=EqGF then
print(`EqGF(S,X,x): the algebraic equation satisfied by the weight-enumerator generating function for the sequence enumerating generalized Dyck paths with set of steps S keeping track of the individual steps. Try:`):
print(`EqGF({1,-1},X,x);`):


elif nops([args])=1 and op(1,[args])=EqGFe then
print(`EqGFe(S,X,x): the algebraic equation satisfied by the weight-enumerator generating function for the sequence enumerating generalized LOSING Dyck paths with set of steps S keeping track of the individual steps. Try:`):
print(`EqGFe({1,-1},X,x);`):


elif nops([args])=1 and op(1,[args])=EqGFt then
print(`EqGFt(S,X,t): the algebraic equation satisfied by X(t), the ordinary generating function  of the sequence enumerating generalized Dyck paths with a set of steps S. Try:`):
print(`EqGFt({1,-1},X,t);`):


elif nops([args])=1 and op(1,[args])=EqGFtS then
print(`EqGFtS(S,X,t): the algebraic equation satisfied by X(t), the ordinary generating function  of the sequence enumerating STRICT (i.e. that never touch the x-axis except at the endpoints), generalized Dyck paths with a set of steps S. Try:`):
print(`EqGFtS({1,-1},X,t);`):

elif nops([args])=1 and op(1,[args])=EqP then
print(`EqP(P,X,t): Given a loaded die P in terms of a list of pairs, outputs the algebraic equation for the duration until being not-in-debt for the last time.`):
print(`to multiply the solution. Try:`):
print(`EqP([[1,1/3],[-1,2/3]],X,t);`):

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

elif nops([args])=1 and op(1,[args])=GFstx then
print(`GFstx(S,x,t,K): The first K terms of the weighted generating functio of the number of generalized Dyck paths with set of steps S with weight of steps x[s]. `):
print(`For checking puroposes only:Try:`):
print(`GFstx({1,0,-1},x,t,6);`):

elif nops([args])=1 and op(1,[args])=GFpr then
print(`GFpr(st,P,t,K): The first K terms of the weighted generating functio of the number of LOSING generalized Dyck paths with set of steps S with probability distribution P. Starting with capital st. Try:`):
print(`GFpr(0,[[1,1/4],[0,1/4],[-1,1/2]],t,6);`):

elif nops([args])=1 and op(1,[args])=IsStrict then
print(`IsStrict(w): Is the walk w that starts at [0,a], stricty above the x-axis except (possibly) the endpoints`):
print(`Try: `):
print(`IsStrict([1,1,-1,-1],0); `):
print(`IsStrict([1,1,-1,-1,1,-1],0);`):


elif nops([args])=1 and op(1,[args])=LD then
print(`LD(P): inputs a list of pairs [face,prob] and outputs face with that prob. Try: `):
print(`LD([[1,1/2],[2,1/3],[-7,1/6]]);`):

elif nops([args])=1 and op(1,[args])=MakeEqF then
print(`MakeEqF(f,g,x,a,b,S): Given symbols f and g  and non-neg. integers a and b`):
print(`outputs an equation for f[a,b]  in terms of other f[a',b'] and g[a',b']. `):
print(`It also returns the quantities used.Try:`):
print(`MakeEqF(f,g,x,0,0,{1,-1});`):


elif nops([args])=1 and op(1,[args])=MakeEqG then
print(`MakeEqG(f,g,x,a,b,S): Given symbols f and g and x a set of integers S, and non-neg. integers a and b`):
print(`outputs an equation for g[a,b] in terms of f[a',b'] and other g[a',b']. `):
print(`It also returns the quantities used.Try:`):
print(`MakeEqG(f,g,x,0,0,{1,2,-1,-2}); `):


elif nops([args])=1 and op(1,[args])=MakeSys then
print(`MakeSys(f,g,x,S): Inputs symbols f, g, and x, and a set of integers S, outouts the`):
print(`system of equations, including f[0,0] needed to find it, along with all the other`):
print(`participating variables. Try:`):
print(`MakeSys(f,g,x,{1,-1});`):



elif nops([args])=1 and op(1,[args])=MakeSysE then
print(`MakeSysE(f,g,x,S,Z): Inputs symbols f, g, and x, and a set of integers S, outputs the`):
print(`system of equations, including a special symbol Z, for  the "gambling histories that winded up losing"`):
print(`along with all the other`):
print(`participating variables. Try:`):
print(`MakeSysE(f,g,x,{1,-1},Z);`):

elif nops([args])=1 and op(1,[args])=MakeSysT then
print(`MakeSysT(f,g,t,S): Inputs symbols f, g, and t, and a set of integers S, outputs the`):
print(`system of equations, including f[0,0] needed to find it, along with all the other`):
print(`participating variables, only keeping track of the length, but not the individual steps.  Try:`):
print(`MakeSysT(f,g,t,{1,-1});`):

elif nops([args])=1 and op(1,[args])=MakeSysT1 then
print(`MakeSysT1(f,g,t,S,st,X): converts MakeSysT(f,g,t,S) to a transoformation format, [LIST,LIST], phrased in terms of an indexed variable X. Try:`):
print(`MakeSysT1(f,g,t,{0,1,-1},f[0,0],X);`):


elif nops([args])=1 and op(1,[args])=ModeP then
print(`ModeP(P,q): inputs a polynomial P in q and outputs the pair [degree,mode] where degree is its degree and mode is the power with the largest coefficient`):
print(`assuming that it is unimodal, it returns FAIL otherwise. Try:`):
print(`ModeP(1+3*q+2*q^2,q);`):
print(`ModeP(qSeqS({1,0,-1},q,10)[11],q);`):

elif nops([args])=1 and op(1,[args])=NuWabn then
print(`NuWabn(a,b,n,S):  the NUMBER of walks from [0,a] to [n,b] using steps [1,s] (s in S) and staying weakly above the x-axis. Try:`):
print(`NuWabn(0,0,5,{0,1,-1});`):


elif nops([args])=1 and op(1,[args])=OneTrip then
print(`OneTrip(P,st): inputs a list of pairs [amount,prob] where the prob. that in one shot you get amount is prob. keeps going until you are in the red, starting with capical st`):
print(`The expected gain of every step must be negative. It outputs the number of steps until you are broke. Try:`):
print(`OneTrip([[1,1/3],[-1,2/3]],0);`):

elif nops([args])=1 and op(1,[args])=Paper1 then
print(`Paper1(S,MaxC): An article about the enumerating sequence of words in the alphabet S that add up to 0 and whose partial sums are never neagtive. It also outputs a linear recurrence  equation satisfied by the`):
print(`coefficients of complexity (ORDER+DEGREE)<=MaxC it one exits.  Try:`):
print(`Paper1({1,0,-1},10):`):


elif nops([args])=1 and op(1,[args])=Paper1A then
print(`Paper1A(S): An article about the enumerating sequence of words in the alphabet S that add up to 0 and the sum of the areas . Try:`):
print(`Paper1A({1,0,-1}):`):


elif nops([args])=1 and op(1,[args])=Paper11A then
print(`Paper11A(S,c): Similar to Paper1A(S), but as numbered theorem c. Try`):
print(`Paper11A({1,0,-1},3):`):



elif nops([args])=1 and op(1,[args])=Paper1p then
print(`Paper1p(P,K1,K2,K3): An article about the algebraic equation satisfied by the probability generating function of  life in a Casino where `):
print(`you are given a finite probability distribution P expressed as  a list of pairs [a,p[a]], where a is an integer and p[a] is a rational number`):
print(`indicating that the probability of getting a in one round is p[a]`):
print(` It also gives you either the exact (if possible), or conditioned upon exiting after L steps`):
print(`and does some confirming simulations with K3 tries. Try:`):
print(`Paper1p([[1,1/3],[0,1/6],[-1,1/2]],4,1000,10000);`):
print(`Paper1p([[2,1/6],[1,1/6],[0,1/6],[-1,1/20],[-2,9/20]],4,200,1000);`):

elif nops([args])=1 and op(1,[args])=Paper1pS then
print(`Paper1pS(k): An article about the algebraic equation satisfied by the SYMBOLIC probability generating function of  life in a Casino where the`):
print(`probability of i (between -k and k) is p[i], and of course they all add-up to 1. Try: `):
print(`Paper1pS(1):`):


elif nops([args])=1 and op(1,[args])=Paper11 then
print(`Paper11(S,MaxC,c): Like Paper1(S,MaxC) but with the theorem named Theorem c, Try:`):
print(`Paper11({1,0,-1},10,3):`):

elif nops([args])=1 and op(1,[args])=Paper2A then
print(`Paper2A(S,MaxC): An article about the sequence giving the number, and sum of areas, of generalized Dyck path with set of steps S and tries to find recurrences`):
print(`for them if their complexity is <=MaxC. Otherwise it just gives the beginning of these sequences. Try:`):
print(`Paper2A({-1,0,1},10);`):


elif nops([args])=1 and op(1,[args])=Paper2pS then
print(`Paper2pS(k,K): An article about the algebraic equation satisfied by the probability generating function of life in a Casino where the`):
print(`probability of i (between  -1 and k) is p[i], and of course they all add-up to 1, and the expected gain is negative, it also gives you `):
print(`explicit expressions for the expectation, standard-deviation, and the 3rd through the k-th moment in termss of the p[i]. try`):
print(`Paper2pS(1,4):`):

elif nops([args])=1 and op(1,[args])=PWabn then
print(`PWabn(a,b,n,P): The PROBABILITY of all the walks from [0,a] to [n,b] using steps [1,s] (s in S) and staying weakly about the x-axis. Try:`):
print(`PWabn(0,0,5,[[1,1/4],[0,1/4],[-1,1/2]]); `):


elif nops([args])=1 and op(1,[args])=PWE then
print(`PWE(st,n,P): The PROBABILITY of a walk with prob. distribution P ending after n steps, starting with capital st. Try:`):
print(`PWE(0,5,[[1,1/4],[0,1/4],[-1,1/2]]); `):

elif nops([args])=1 and op(1,[args])=PosW then
print(`PosW(S,n): Given a set of integers (both positive and negative) and a non-negative integer  n outputs all the`):
print(`words in S such that the partial sums are >=0 and whose sum is 0. Try:`):
print(`PosW({0,1,-1},5);`):

elif nops([args])=1 and op(1,[args])=PosW1 then
print(`PosW1(S,n,k): Given a set of integers (both positive and negative) and non-negative integer  n and k outputs all the`):
print(`words in S such that the partial sums are >=0 and whose sum is k. Try:`):
print(`PosW1({0,1,-1},5,1);`):


elif nops([args])=1 and op(1,[args])=qEqGFt then
print(`qEqGFt(S,X,t): Inputs a set of positive and negative integers S, outputs an algebraic equation for X=X(t), the`):
print(`generating function for the sum-of-the-areas under ALL generalized Dyck paths with set of steps {[1,s], s in S}. For example, for the the`):
print(`classical Dyck case do:`):
print(`qEqGFt({1,-1},X,t);`):
print(`For Motzkin walks, enter:`):
print(`qEqGFt({1,0,-1},X,t);`):


elif nops([args])=1 and op(1,[args])=qEqGFtS then
print(`qEqGFtS(S,X,t): Inputs a set of positive and negative integers S, outputs an algebraic equation for X=X(t), the`):
print(`generating function for the sum-of-the-areas under ALL STRICT (i.e. that do not touch the x-axis except for the endpoints) generalized Dyck paths with set of steps {[1,s], s in S}. For example, for the the`):
print(`classical Dyck case do:`):
print(`qEqGFt({1,-1},X,t);`):
print(`For Motzkin walks, enter:`):
print(`qEqGFtS({1,0,-1},X,t);`):


elif nops([args])=1 and op(1,[args])=qMakeEqFt then
print(`qMakeEqFt(f, g, t, q,a, b, S): Sets up a functional equation for the quantity f[a,b](t,q), the set of walks using the set of steps {[1,s], s in S}, that start at [0,a] and end where y=b`):
print(`and that are WEAKLY above the x-axis. Try:`):
print(`qMakeEqFt(f,g,t,q,0,0,{1,0,-1});`):


elif nops([args])=1 and op(1,[args])=qMakeEqGt then
print(`qMakeEqGt(f, g, t, q,a, b, S): Sets up a functional equation for the quantity g[a,b](t,q), the set of walks using the set of steps {[1,s], s in S}, that start at [0,a] and end where y=b`):
print(`and that are STRICTLY above the x-axis (expcept, possibly at the endpoints). Try:`):
print(`qMakeEqGt(f,g,t,q,0,0,{1,0,-1});`):

elif nops([args])=1 and op(1,[args])=qMakeSyst then
print(`qMakeSyst(f,g,t,q,S): inputs symbols f, g, for functions, and t,q, for variables, and a set of (positive and negative integers).`):
print(`Output: the  a system of`):
print(`functional equations for the quantities f[a,b](t,q) and g[a,b](t,q), where `):
print(`f[a,b](t,q) (respectively g[a,b](t,q)) is the weight-enumerator according to the weight`):
print(`t^length(w)*q^area(w) of all paths using the steps {[1,s], s in S} starting at [0,a] and ending at [n,b] and WEAKLY (resp. STRICTLY, except possibly at the endpoints) above the x-axis.`):
print(`It also outputs the set of functions f[a,b](t), f[a,b](q*t), g[a,b](t), for various a and b, that feature in the equations. It is the minimal system that enables us to deterine our goal-in-life, f[0,0](t,q). The dependence on q`):
print(`is implied, so, all quantities also depend on the variable q. `):
print(``):
print(`For example, the system for Motzkin paths is gotten by typing:`):
print(`qMakeSyst(f,g,t,q,{1,0,-1});`):

elif nops([args])=1 and op(1,[args])=qSeqS then
print(`qSeqS(S,q,K): The first K terms, starting with n=0 of the weight-enumerator, according to q^area,  of generalized Dyck path with set steps S .`):
print(`Try:`):
print(`qSeqS({1,0,-1},q,10);`):

elif nops([args])=1 and op(1,[args])=SeqFromEq then
print(`SeqFromEq(Eq,X,t,K): Given an equation Eq satisfied by X=X(t) finds the first K coefficients. Try:`):
print(`SeqFromEq(t*X^2+1-X,X,t,10);`):

elif nops([args])=1 and op(1,[args])=SeqFromSysEq then
print(`SeqFromSysEq(Eqs,t,X,K): Given a system of  equations in terms of X[1],X[2],... and a variable T`):
print(`finds the first K terms of X[1]. Try:`):
print(`gu:=MakeSysT1(f,g,t,{1,-1,0},f[0,0],X): SeqFromSysEq(gu,t,X,10);`):


elif nops([args])=1 and op(1,[args])=SeqS then
print(`SeqS(S,K): The first K+1 terms, starting at n=0 of the sequence enumerating sequences of length n in the alphabet S (of integers) with the property that the sum is 0 and the partial sums are never negative`):
print(`In other words generalized Dyck words except that instead of the alphabet {1,-1} you have the alphabet S. DONE via numerical dynamical programming. For example to get the`):
print(`first 31 Motzkin numbers, starting at n=0, type:`):
print(`SeqS({1,0,-1},30);`):

elif nops([args])=1 and op(1,[args])=SeqSe then
print(`SeqSe(S,K): Same as SeqS(S,K) but using the equation satisfied by the generating function obtained via SeqFromEq(EqGFt(S,X,t),X,t,K). Try`):
print(`SeqSe({1,0,-1},30);`):


elif nops([args])=1 and op(1,[args])=SeqSee then
print(`SeqSee(S,K): Same as SeqS(S,K) and SeqSe(S,K) but using the SYSTEM of equation satisfied by the generating function obtained via MakeSysT1. Try`):
print(`SeqSee({1,0,-1},30);`):

elif nops([args])=1 and op(1,[args])=SimuGF then
print(`SimuGF(P,x,K,st): The empirical generating function of the r.v. "number of rolls until you are in the red" with loaded die P, using K sessions, starting with capital st. Try:`):
print(`SimuGF([[1,1/3],[0,1/3],[-2,1/3]],x,1000,0);`):


elif nops([args])=1 and op(1,[args])=SimuStat then
print(`SimuStat(P,K1,N,st): Given a loaded losing rambling die P, gambles K1 times and then outputs the average, standard deviation, and the scaled moments up to the N-th.  The starting capital is st.`):
print(`It also outputs the max. longevity in these K1 runs`):
print(`Try:`):
print(`SimuStat([[1,1/3],[0,1/3],[-2,1/3]],1000,4,0);`):

elif nops([args])=1 and op(1,[args])=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 nops([args])=1 and op(1,[args])=WabnSdirect then
print(`WabnSdirect(a,b,n,S): Same as WabnS(a,b,n,S) but done diretly. For checking purpose only. Try:`):
print(`WabnSdirect(0,0,5,{0,1,-1}):`):

elif nops([args])=1 and op(1,[args])=Wabn then
print(`Wabn(a,b,n,S): The SET of walks from [0,a] to [n,b] using steps [1,s] (s in S) and staying weakly above the x-axis. Try:`):
print(`Wabn(0,0,5,{0,1,-1});`):


elif nops([args])=1 and op(1,[args])=WabnS then
print(`WabnS(a,b,n,S): all the walks from [0,a] to [n,b] using steps [1,s] (s in S) and staying STIRCTLY above the x-axis, excpet possibly the end points. Try:`):
print(`WabnS(0,0,5,{0,1,-1});`):


elif nops([args])=1 and op(1,[args])=Wt then
print(`Wt(w,x): The weight of the walk  using x[s] for a step s. Try:`):
print(`Wt([1,0,2,-1],x);`):


print(``):
else
 print(`There is no such thing as`, args):

fi:


end:



###From FindRec
ezraF:=proc()
if args=NULL then

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

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

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

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


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


fi:
 
end:


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

ope:=subs(var1,ope):

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

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

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

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


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

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

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


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

ope:=subs(var1,ope):

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

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

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

end:
 


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


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

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:

 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):

kak:=denom(normal(ope)):

if member(0,{seq(subs(n=i,kak),i=1..100)}) then
  RETURN(FAIL):
fi:

if ope<>FAIL  then

if  SeqFromRec(ope,n,N,[op(1..degree(ope,N),f)],nops(f))=f then
 RETURN(ope):
else
 RETURN(FAIL):
fi:

fi:
 od:
od:
FAIL:

end:





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

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

gu:=Ini:

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

gu:

end:
 
#end Findrec


with(linalg):

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


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


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

A1:=matrix(D1,D1):

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

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

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

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

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

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

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

true:

end:
###End From FindRec

#PosW1(S,n,k): Given a set of integers (both positive and negative) and non-negative integer  n and k outputs all the
#words in S such that the partial sums are >=0 and whose sum is k. Try:
#PosW1({0,1,-1},5,1);
PosW1:=proc(S,n,k) local gu,lu,lu1,s:
option remember:
if n<0 or k<0 then
 RETURN({}):
fi:

if n=0 then
 if k=0 then
  RETURN({[]}):
 else
 RETURN({}):
fi:
fi:

gu:={}:

for s in S do
 lu:=PosW1(S,n-1,k-s):
 gu:=gu union {seq([op(lu1),s], lu1 in lu)}:
od:

gu:
end:



#PosW(S,n): Given a set of integers (both positive and negative) and non-negative integer  n outputs the set of
#words in S such that the partial sums are >=0 and whose sum is 0. Try:
#PosW({0,1,-1},5);
PosW:=proc(S,n):PosW1(S,n,0) :end:


#Wabn(a,b,n,S): all the walks from [0,a] to [n,b] using steps [1,s] (s in S) and staying weakly about the x-axis. Try:
#Wabn(0,0,5,{0,1,-1}):
Wabn:=proc(a,b,n,S) local gu,lu,lu1,s:
option remember:
if a<0 or b<0 then
RETURN({}):
fi:

if n=0 then
 if a=b then
  RETURN({[]}):
else
 RETURN({}):
fi:
fi:

gu:={}:


for s in S do
 lu:=Wabn(a,b-s,n-1,S):
 gu:=gu union {seq([op(lu1),s], lu1 in lu)}:
od:


gu:

end:




#NuWabn(a,b,n,S): The NUMEBR of all the walks from [0,a] to [n,b] using steps [1,s] (s in S) and staying weakly about the x-axis. Try:
#NuWabn(0,0,5,{0,1,-1}):
NuWabn:=proc(a,b,n,S) local gu,lu,lu1,s:
option remember:
if a<0 or b<0 then
RETURN(0):
fi:

if n=0 then
 if a=b then
  RETURN(1):
else
 RETURN(0):
fi:
fi:

add(NuWabn(a,b-s,n-1,S),s in S):

end:



#PWabn(a,b,n,P): The PROBABILITY of all the walks from [0,a] to [n,b] using steps [1,s] (s in S) and staying weakly about the x-axis. Try:
#PWabn(0,0,5,{[[1,1/4],[0,1/4],[-1,1/2]}):
PWabn:=proc(a,b,n,P) local gu,lu,lu1,s,i:
option remember:
if a<0 or b<0 then
RETURN(0):
fi:

if n=0 then
 if a=b then
  RETURN(1):
else
 RETURN(0):
fi:
fi:

add(P[i][2]*PWabn(a,b-P[i][1],n-1,P),i=1..nops(P)):

end:



#PWE(st,n,P): The PROBABILITY of a walk with prob. distribution P,ending after n steps, starting with capital st.Try:
#PWE(0,5,[[1,1/4],[0,1/4],[-1,1/2]]):
PWE:=proc(st,n,P) local gu,i,s,a:
gu:=0:
for i from 1 to nops(P) do
s:=P[i][1]:
if s<0 then
for a from  0 to -s-1 do
gu:=gu+P[i][2]*PWabn(st,a,n-1,P):
od:
fi:
od:
gu:
end:



#IsStrict(w): Is the walk w that starts at [0,a], stricty above the x-axis except (possibly) the endpoints
#Try:
#IsStrict([1,1,-1,-1],0);
#IsStrict([1,1,-1,-1,1,-1],0);

IsStrict:=proc(w,a) local i,s:
s:=a:
for i from 1 to nops(w)-1 do
s:=s+w[i]:
if s=0 then
 RETURN(false):
fi:
od:
true:
end:





#WabnSdirect(a,b,n,S): all the Strict walks from [0,a] to [n,b] using steps [1,s] (s in S) and staying stricly above the x-axis. Try:
#WabnSdirect(0,0,5,{0,1,-1}):
WabnSdirect:=proc(a,b,n,S) local lu,lu1,gu:
lu:=Wabn(a,b,n,S):
gu:={}:

for lu1 in lu do
 if IsStrict(lu1,a) then
  gu:=gu union {lu1}:
 fi:
od:
gu:
end:






#WabnS1(a,b,n,S): all the STRICT walks from [0,a] to [n,b] using steps [1,s] (s in S) and staying strictly above the x-axis. Try:
#WabnS1(0,0,5,{0,1,-1}):
WabnS1:=proc(a,b,n,S) local gu,lu,lu1,s:
option remember:
if n=0 then
 if a=b then
  RETURN({[]}):
else
 RETURN({}):
fi:
fi:

gu:={}:


if a<0 or b<=0 then
RETURN({}):
fi:

for s in S do
 lu:=WabnS1(a,b-s,n-1,S):
 gu:=gu union {seq([op(lu1),s], lu1 in lu)}:
od:


gu:

end:




#WabnS(a,b,n,S): all the STRICT walks from [0,a] to [n,b] using steps [1,s] (s in S) and staying strictly above the x-axis expcet possibly the end points. Try:
#WabnS(0,0,5,{0,1,-1}):
WabnS:=proc(a,b,n,S) local gu,lu,lu1,s:
gu:={}:

if n=0 then
RETURN({}):
fi:

for s in S do

 lu:=WabnS1(a,b-s,n-1,S):
 gu:=gu union {seq([op(lu1),s], lu1 in lu)}:

od:
gu:
end:


#DecoW(w,a): Given a walk that starts at [0,a] outputs [w1,w2] where w1 is the  largest prefix that is strictly above the x-axis, followed by the rest so w=w1w2. Try
#DecoW([1,-1,1,-1],0);
DecoW:=proc(w,a) local s,i:

s:=a+w[1]:


for i from 2 to nops(w) while s>0 do
 s:=s+w[i]:
od:
[[op(1..i-1,w)],[op(i..nops(w),w)]]:

end:





#MakeEqF(f,g,x,a,b,S): Given symbols f and g  and non-neg. integers a and b
#outputs what f[a,b] is equal to in terms of outer f[a1,b1] and g[a1,b1]. 
#It also returns the quantities used
#Try:
#MakeEqF(f,g,x,0,0,S):
MakeEqF:=proc(f,g,x,a,b,S):

if a>0 and b>0 then
RETURN(f[a,b]=f[a-1,b-1]+g[a,0]*f[0,b], {f[a-1,b-1],g[a,0],f[0,b]}):
fi:

if a=0 and b>0 then
RETURN(f[a,b]=f[0,0]*g[0,b], {f[0,0],g[0,b]}):
fi:

if a>0 and b=0 then
 RETURN(f[a,b]=g[a,0]*f[0,0],{g[a,0],f[0,0]}):
fi:

if a=0 and b=0 then
 if not member(0,S) then
    RETURN(f[0,0]=1+f[0,0]*g[0,0],{f[0,0],g[0,0]}):
  else
   RETURN(f[0,0]=1+x[0]*f[0,0]+f[0,0]*g[0,0],{f[0,0],g[0,0]}):
 fi:
fi:

end:



#MakeEqG(f,g,x,a,b,S): Given symbols f and g and x a set of integers S, and non-neg. integers a and b
#outputs what g[a,b] is equal to in terms of f[a',b'] and other g[a',b']. Try
#MakeEqG(f,g,x,0,0,{1,2,-1,-2}):
MakeEqG:=proc(f,g,x,a,b,S) local gu,s,P,N,s1,s2,lu:
P:={}:
N:={}:



for s in S do
 if s>0 then
 P:=P union {s}:
elif s<0 then
 N:=N union {s}:
fi:
od:


if a>0 and b>0 then
print(`Something is wrong`):
 RETURN(FAIL):
fi:

gu:=0:
lu:={}:
if a=0 and b>0 then
for s  in P do
 lu:=lu union {f[a+s-1,b-1]}:
 gu:=gu+x[s]*f[a+s-1,b-1]:
od:
RETURN(g[a,b]=gu,lu):
fi:


if a>0 and b=0 then
for s  in N do
 lu:=lu union {f[a-1,b-s-1]}:
gu:=gu+x[s]*f[a-1,b-s-1]:
od:
RETURN(g[a,b]=gu,lu):
fi:

if a=0 and b=0 then

for s1 in P do
 for s2 in N do
lu:=lu union {f[a+s1-1,b-s2-1]}:
  gu:=gu+x[s1]*x[s2]*f[a+s1-1,b-s2-1]:
 od:
od:

RETURN(g[0,0]=gu,lu):
fi:

end:


#MakeSys(f,g,x,S): Inputs symbols f, g, and x, and a set of integers S, outouts the
#system of equations, including f[0,0] needed to find it, along with all the other
#participating variables. Try:
#MakeSys(f,g,x,{1,-1});
MakeSys:=proc(f,g,x,S) local eq,StillToDo,AlreadyDone,guy,ku:
eq:={}:
StillToDo:={f[0,0]}:

AlreadyDone:={}:

while StillToDo<>{} do
guy:=StillToDo[1]:

if op(0,guy)=f then
 ku:=MakeEqF(f,g,x,op(1,guy),op(2,guy),S):
 AlreadyDone:=AlreadyDone union {guy}:
 StillToDo:=StillToDo union ku[2] minus AlreadyDone:
 eq:=eq union {ku[1]}:

 elif op(0,guy)=g then
 ku:=MakeEqG(f,g,x,op(1,guy),op(2,guy),S):
 AlreadyDone:=AlreadyDone union {guy}:
 StillToDo:=StillToDo union ku[2] minus AlreadyDone:
 eq:=eq union {ku[1]}:
else
 print(`Something is wrong`):
 RETURN(FAIL):
fi:
od:
eq,AlreadyDone:

end:




#MakeSysT(f,g,t,S): Inputs symbols f, g, and t, and a set of integers S, outouts the
#system of equations, including f[0,0] needed to find it, along with all the other
#participating variables, only keeping track of the length, but not the individual steps.  Try:
#MakeSysT(f,g,t,{1,-1});
MakeSysT:=proc(f,g,t,S) local gu,x,i,eq,var,s1:
gu:=MakeSys(f,g,x,S):
eq:=gu[1]:
var:=gu[2]:
eq:=subs({seq(x[s1]=t,s1 in S)},eq):
eq,gu[2]:
end:




#EqGF(S,X,x): the algebraic equation satisfied by the generating function for the sequence enumerating generalized Dyck paths with set of steps S. Try:
#EqGF({1,-1},X,x);
EqGF:=proc(S,X,x) local gu,f,g,eq,var,eq1:
gu:=MakeSys(f,g,x,S):
eq:=gu[1]:
var:=gu[2] minus {f[0,0]}:

eq:={seq(op(1,eq1)-op(2,eq1),eq1 in eq)}:
subs(f[0,0]=X,Groebner[Basis](eq,plex(op(var),f[0,0]))[1]):

end:



#EqGFt(S,X,t): the algebraic equation satisfied by the generating function for the sequence enumerating generalized Dyck paths with set of steps S of length n not keeping track of the original steps. Try:
#EqGFt({1,-1},X,t);
EqGFt:=proc(S,X,t) local gu,x,f,g,eq,var,s1,eq1,lu,Neq,i:
option remember:
gu:=MakeSys(f,g,x,S):
eq:=gu[1]:
var:=gu[2] minus {f[0,0]}:

eq:=subs({seq(x[s1]=t,s1 in S)},eq):


eq:={seq(op(1,eq1)-op(2,eq1),eq1 in eq)}:

eq:=subs(f[0,0]=X,Groebner[Basis](eq,plex(op(var),f[0,0]))[1]):
eq:=factor(eq):

if not type(eq,`*`) then
eq:=add(factor(coeff(eq,X,i))*X^i,i=0..degree(eq,X)):
 RETURN(eq):
fi:

lu:=add(NuWabn(0,0,i,S)*t^i,i=0..31):

for i from 1 to nops(eq) do
Neq:=op(i,eq):
if {seq(coeff(subs(X=lu,Neq),t,i),i=0..30)}={0} then
Neq:=add(factor(coeff(Neq,X,i))*X^i,i=0..degree(Neq,X)):
 RETURN(Neq):
fi:
od:



end:



#Wt(w,x): The weight of the walk  using x[s] for a step s. Try:
#Wt([1,0,2,-1],x);
Wt:=proc(w,x) local i:
 mul(x[w[i]],i=1..nops(w)):
end:

#GFstx(S,x,t,K): The first K terms of the weighted generating functio of the number of generalized Dyck paths with set of steps S with weight of steps x[s]. 
#For checking puroposes only:Try:
#GFstx({1,0,-1},x,t,6);
GFstx:=proc(S,x,t,K) local gu,lu,i,lu1:
gu:=0:
for i from 0 to K do
lu:=Wabn(0,0,i,S):
lu:=add(Wt(lu1,x),lu1 in lu):
gu:=gu+t^i*lu:
od:
gu:

end:


#GFprold(P,t,K): The first K terms of the weighted generating functio of the number of LOSING generalized Dyck paths with set of steps S with probability distribution P
#For checking puroposes only:Try:
#GFprold([[1,1/4],[0,1/4],[-1,1/2]],x,t,6);
GFprPold:=proc(P,t,K) local gu,lu,i,lu1,S,s,kv,kv1,a,x:
S:={seq(P[i][1],i=1..nops(P))}:
gu:=0:
for i from 0 to K do
lu:=0:
for s in S do
 if s<0 then
  for a from 0 to -s-1 do
   kv:=Wabn(0,a,i,S):
 lu:=lu+add(Wt(kv1,x),kv1 in kv)*x[s]:
od:
 fi:
od:
gu:=gu+subs({seq(x[P[i][1]]=P[i][2]*t,i=1..nops(P))},lu):
od:
gu:

end:




#GFpr(st,P,t,K): The first K terms of the weighted generating functio of the number of LOSING generalized Dyck paths with set of steps S with probability distribution P
#starting with capital st. Try:
#GFpr(0,[[1,1/4],[0,1/4],[-1,1/2]],x,t,6);
GFpr:=proc(st,P,t,K) local i:
add(PWE(st,i,P)*t^i,i=1..K+1):
end:



#CheckEqGF(S,K): Checks procedure EqGF(S,X,x) up to K terms. Try:
#CheckEqGF({1,0,-1},10);
CheckEqGF:=proc(S,K) local X,t,lu,gu,i,x,s:
gu:=EqGF(S,X,x) :
gu:=subs({seq(x[s]=x[s]*t,s in S)},gu):
lu:=GFstx(S,x,t,K+1):
gu:=expand(subs(X=lu,gu)):

evalb([seq(expand(coeff(gu,t,i)),i=0..K)]=[0$(K+1)]):
end:



#CheckEqGFt(S,K): Checks procedure EqGFt(S,X,t) up to K terms. Try:
#CheckEqGFt({1,0,-1},10);
CheckEqGFt:=proc(S,K) local X,t,lu,gu,i:
gu:=EqGFt(S,X,t) :
lu:=add(NuWabn(0,0,i,S)*t^i,i=0..K):
gu:=expand(subs(X=lu,gu)):

evalb([seq(expand(coeff(gu,t,i)),i=0..K)]=[0$(K+1)]):
end:



#CheckEqP(P,K): Checks procedure EqP(P,X,t) up to K terms. Try:
#CheckEqP({1,0,-1},10);
CheckEqP:=proc(P,K) local X,t,lu,gu,i:
gu:=EqP(P,X,t) :
lu:=GFpr(0,P,t,K):
gu:=expand(subs(X=lu,gu)):
evalb({seq(coeff(gu,t,i),i=1..K)}={0}):
end:



#SeqFromEq(Eq,X,t,K): Given an equation Eq satisfied by X=X(t) finds the first K coefficients. Try:
#SeqFromEq(t*X^2+1-X,X,t,10);
SeqFromEq:=proc(Eq,X,t,K) local gu,lu,f1,f2,f3,i,i1:


lu:=coeff(coeff(Eq,X,0),t,0):

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

gu:=normal(Eq/lu+X):


f1:=1:
f2:=expand(subs(X=f1,gu)):
f2:=add(coeff(f2,t,i1)*t^i1,i1=0..K):


for i from 1 to K while f1<>f2 do
 f3:=expand(subs(X=f2,gu)):
 f3:=add(coeff(f3,t,i1)*t^i1,i1=0..K):
f1:=f2:
f2:=f3:
od:

if i=K+1 and f1<>f2 then
 FAIL:
else
[seq(coeff(f2,t,i1),i1=0..K)]:
fi:

end:

#SeqS(S,K): The first K+1 terms, starting at n=0 of the sequence enumerating sequences of length n in the alphabet S (of integers) with the property that the sum is 0 and the partial sums are never negative
#In other words generalized Dyck words except that instead of the alphabet {1,-1} you have the alphabet S. DONE via numerical dynamical programming. For example to get the
#first 31 Motzkin numbers, starting at n=0, type:
#SeqS({1,0,-1},30);
SeqS:=proc(S,K) local i:
[seq(NuWabn(0,0,i,S),i=0..K)]:
end:


#SeqSe(S,K): Same as SeqS(S,K) but using the equation satisfied by the generating function obtained via SeqFromEq(EqGFt(S,X,t),X,t,K). Try
#SeqSe({1,0,-1},30);
SeqSe:=proc(S,K) local X,t:
SeqFromEq(EqGFt(S,X,t),X,t,K):
end:


#Paper1(S,MaxC): An article about the enumerating sequence of words in the alphabet S that add up to 0 and whose partial sums are never neagtive. Try:
#Paper1({1,0,-1}):

Paper1:=proc(S,MaxC) local X,t,eq,L,n,ope,a,N,i,L1,gadol,L2,t0:
t0:=time():
print(``):
print(`Enumerating generalized Dyck words in the alphabet`, S):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
eq:=EqGFt(S,X,t):
print(``):
print(`Theorem:  Let a(n) be the number of words of length n in the alphabet`, S, ` that sum-up to 0 and whose partial sums are never negative, in other words generalized Dyck words with alphabet`, S, `rather than {1,-1} `):
print(``):
print(`Let X(t) be the ordinary generating function of that sequence, in other words`):
print(``):
print(X(t)=Sum(a(n)*t^n,n=0..infinity)):
print(``):
print(`X(t) satisfies the algebraic equation `):
print(``):
print(subs(X=X(t),eq)=0):
print(``):
print(`and in Maple notation `):
print(``):
lprint(subs(X=X(t),eq)=0):
print(``):
print(`For the sake of the OEIS here are the first 40 terms, starting with n=0`):
print(``):
lprint(SeqS(S,40)):
print(``):
L:=SeqS(S,max(seq((2+i)*(1+MaxC-i),i=0..MaxC))+11):
L:=[op(2..nops(L),L)]:

ope:=Findrec(L,n,N,MaxC):



if ope<>FAIL then

gadol:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),L)],1000)[1000]:

ope:=expand(1-subs(n=n-degree(ope,N),ope)/N^degree(ope,N)):

print(`Furthermore a(n) satisfies the following linear recurrence equation with polynomial coefficients `):
print(``):
print(a(n)=add(factor(coeff(ope,N,-i))*a(n-i),i=1..-ldegree(ope,N))):
print(``):
print(`subject to the initial conditions`):
print(``):
print(seq(a(i)=L[i],i=1..-ldegree(ope,N))):
print(``):
print(`and in Maple notation `):
print(``):
lprint(a(n)=add(factor(coeff(ope,N,-i))*a(n-i),i=1..-ldegree(ope,N))):
print(``):
print(``):
lprint(seq(a(i)=L[i],i=1..-ldegree(ope,N))):
print(``):
print(`Just for fun, using this recurrence we get that`):
print(``):
print(a(1000)=gadol):
print(``):
fi:

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

end:





#Paper11(S,MaxC,c): An article about the enumerating sequence of words in the alphabet S that add up to 0 and whose partial sums are never neagtive. Try:
#Paper11({1,0,-1},10,3):

Paper11:=proc(S,MaxC,c) local X,t,eq,L,n,ope,a,N,i,L1,gadol:
print(``):
print(``):
eq:=EqGFt(S,X,t):
print(``):
print(cat(`Theorem `,c,` :`)):
print(`Let a(n) be the number of words of length n in the alphabet`, S, ` that sum-up to 0 and whose partial sums are never negative, in other words generalized Dyck words with alphabet`, S, `rather than {1,-1} `):
print(``):
print(`Let X(t) be the ordinary generating function of that sequence, in other words`):
print(``):
print(X(t)=Sum(a(n)*t^n,n=0..infinity)):
print(``):
print(`X(t) satisfies the algebraic equation `):
print(``):
print(subs(X=X(t),eq)=0):
print(``):
print(`and in Maple notation `):
print(``):
lprint(subs(X=X(t),eq)=0):
print(``):
print(`For the sake of the OEIS here are the first 40 terms, starting with n=0`):
print(``):
lprint(SeqS(S,40)):
print(``):
L:=SeqS(S,max(seq((2+i)*(1+MaxC-i),i=0..MaxC))+11):
L:=[op(2..nops(L),L)]:
ope:=Findrec(L,n,N,MaxC):
print(``):
if ope<>FAIL then
gadol:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),L)],1000)[1000]:
 L1:=[op(2..100,SeqS(S,100))]:

if L1<>SeqFromRec(ope,n,N,[op(1..degree(ope,N),L1)],nops(L1)) then
print(`Something terrible happened `):
 RETURN(FAIL):
fi:


ope:=expand(1-subs(n=n-degree(ope,N),ope)/N^degree(ope,N)):

print(`Furthermore a(n) satisfies the following linear recurrence equation with polynomial coefficients `):
print(``):
print(a(n)=add(factor(coeff(ope,N,-i))*a(n-i),i=1..-ldegree(ope,N))):
print(``):
print(`subject to the initial conditions`):
print(``):
print(seq(a(i)=L[i],i=1..-ldegree(ope,N))):
print(``):
print(`and in Maple notation `):
print(``):
lprint(a(n)=add(factor(coeff(ope,N,-i))*a(n-i),i=1..-ldegree(ope,N))):
print(``):
print(``):
lprint(seq(a(i)=L[i],i=1..-ldegree(ope,N))):
print(``):
print(`Just for fun, using this recurrence we get that`):
print(``):
print(a(1000)=gadol):
print(``):

fi:



end:


 


#Book(M,MaxC): Inputs a positive integer M and outputs a book about enumearing generalized Dyck path with all possible alphabets consisting of both negative and positive integes and 0. It also
#outputs recurrences if their complexity is <=MaxC. Try:
#Book(2,10);
Book:=proc(M,MaxC) local NEG,POS,c,i,gu1,gu2,gu11,gu21,t0,S:
t0:=time():
print(``):
print(`Enumerating Generalized Dyck paths with alphabets consisting of integers from`, -M, ` to `, M):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
c:=0:
NEG:={seq(-i,i=1..M)}:
POS:={seq(i,i=1..M)}:
gu1:=powerset(NEG) minus {{}}:
gu2:=powerset(POS) minus {{}}:
print(``):
print(`-------------------------------------------------------------`):
print(``):

for gu11 in gu1 do
for gu21 in gu2 do
 S:=gu11 union gu21:

if igcd(op(S))=1 then

c:=c+1:

Paper11(S,MaxC,c):
print(``):
print(`-------------------------------------------------------------`):
print(``):
S:=S union {0}:
c:=c+1:
Paper11(S,MaxC,c):
fi:
od:
od:

print(``):
print(`--------------------------------`):
print(``):
print(`This ends this book that took`, time()-t0, `seconds to create `):
print(``):
end:


#start simulation

#LD(P): inputs a list of pairs [face,prob] and outputs face with that prob. Try
#LD([[1,1/2],[2,1/3],[-7,1/6]]);
LD:=proc(P) local i,d,P1,tot,ra,j:

P1:=[seq(P[i][2],i=1..nops(P))]:

if not (convert(P1,`+`)=1  and min(op(P1))>=0) then
 print(P1, `is not a prob. dist. `):
 RETURN(FAIL):
fi:
 
d:=lcm(seq(denom(P1[i]),i=1..nops(P1))):
P1:=[seq(d*P1[i],i=1..nops(P1))]:
tot:=convert(P1,`+`):
ra:=rand(1..tot)():


for i from 1 to nops(P1) while add(P1[j],j=1..i)<ra do od:
P[i][1]:

end:

#OneTrip(P,st): inputs a list of pairs [amount,prob] where the prob. that in one shot you get amount is prob. keeps going until you are in the red, starting with capical st
#The expected gain of every step must be negative. It outputs the number of steps until you are broke. Try:
#OneTrip([[1,1/3],[-1,2/3]],0);

OneTrip:=proc(P,st) local su,i:

if add(P[i][1]*P[i][2],i=1..nops(P))>=0 then
 print(P, `is not a losing roullete`):
 RETURN(FAIL):
fi:

su:=st:


for i from 1 while su>=0 do
 su:=su+LD(P):
od:
i-1:
end:

#SimuGF(P,x,K,st): The empirical generating function of the r.v. "number of rolls until you are in the red" with loaded die P, using K sessions, starting with capital st. Try:
#SimuGF([[1,1/3],[0,1/3],[-2,1/3]],x,1000,0);
SimuGF:=proc(P,x,K,st) local i:

if not (type(P,list) and add(P[i][2],i=1..nops(P))=1 and min(seq(P[i][2],i=1..nops(P)))>=0 and add(P[i][1]*P[i][2],i=1..nops(P))<0 )then
print(`bad input`):
 RETURN(FAIL):
fi:

add(x^OneTrip(P,st),i=1..K):
end:


#Stat(F,t,K): Given a generating function F of t, for a r.v. finds the statistical information. Try:
#Stat((1+t)^10000,t,4);
Stat:=proc(F,t,K) local gu,mu,sd,k,f:
f:=evalf(F/subs(t=1,F)):
mu:=subs(t=1,diff(f,t)):
gu:=[mu]:
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:
gu:

end:

#SimuStat(P,K1,N): Given a loaded losing rambling die P, gambles K1 times and then outputs the average, standard deviation, and the scaled moments up to the N-th.  Starting with capital st,
#It also outputs the max. longevity in these K1 runs
#Try
#SimuStat[[1,1/3],[0,1/3],[-2,1/3]],1000,4,0);
SimuStat:=proc(P,K1,N,st) local x,i,f:

if not (type(P,list) and add(P[i][2],i=1..nops(P))=1 and min(seq(P[i][2],i=1..nops(P)))>=0 and add(P[i][1]*P[i][2],i=1..nops(P))<0 and st>=0 )then
print(`bad input`):
 RETURN(FAIL):
fi:

f:=SimuGF(P,x,K1,st):
Stat(f,x,N),degree(f,x):
end:




#MakeSysE(f,g,x,S,Z): Inputs symbols f, g, and x, and a set of integers S, outouts the
#system of equations, including a special symbol Z, for  the "gambling histories that winded up losing"
#along with all the other
#participating variables. Try:
#MakeSysE(f,g,x,{1,-1},Z);
MakeSysE:=proc(f,g,x,S,Z) local gu,eq,var,eq1,a,s:

gu:=MakeSys(f,g,x,S):
eq:=gu[1]:
var:=gu[2] union {Z}:

eq1:=0:

for s in S do
 if s<0 then
  for a from 0 to -s-1 do
  eq1:=eq1+f[0,a]*x[s]:
  od:
 fi:
od:

eq:=eq union {Z=eq1}:


eq,var:

end:


#EqGFe(S,X,x): the algebraic equation satisfied by the generating function for the sequence enumerating LOSING generalized Dyck paths with set of steps S. Try:
#EqGFe({1,-1},X,x);
EqGFe:=proc(S,X,x) local gu,f,g,eq,var,eq1:
gu:=MakeSysE(f,g,x,S,X):
eq:=gu[1]:
var:=gu[2] minus {X}:

eq:={seq(op(1,eq1)-op(2,eq1),eq1 in eq)}:
Groebner[Basis](eq,plex(op(var),X))[1]:

end:


#EqP(P,X,t): Given a loaded die P in terms of a list of pairs, outputs the algebraic equation for the duration until being not-in-debt for the last time
#to multiply the solution. Try:
#EqP([[1,/3],[-1,2/3]],X);
EqP:=proc(P,X,t) local x,S,f,i:

S:={seq(P[i][1],i=1..nops(P))}:

f:=EqGFe(S,X,x):
subs({seq(x[P[i][1]]=P[i][2]*t,i=1..nops(P))},f):
end:



###start new stuff




#DerF(F,X,t): inputs a polynomial f in X and t, uses implicit differention to finds an explicit expression for X'(t), in terms of X and t where X stands for X(t)
#Try: DerF(X-1-t*X^2,X,t);
DerF:=proc(F,X,t):
-normal(diff(F,t)/diff(F,X)):
end:


#DerEqP(P,X,t): an expression for X'(t) if X(t) is the prob. gen. function for the Genralized Dyck path process. Try:
#DerEqP([[1,a],[0,b],[-1,1-a-b]],X,t);
DerEqP:=proc(P,X,t) local gu,c1,gu1,d,lu:
gu:=EqP(P,X,t):
d:=degree(gu,X):
c1:=coeff(gu,X,d):
gu1:=normal((c1*X^d-gu)/c1):
lu:=DerF(gu,X,t):
lu:=normal(t*subs(X^d=gu1,lu)):
lu:
end:


#DerEqPk(P,X,t,k): an expression for (t*d/dt)^k X(t) if X(t) is the prob. gen. function for the Genralized Dyck path process. Try:
#DerEqPk([[1,a],[0,b],[-1,1-a-b]],X,t,2);
DerEqPk:=proc(P,X,t,k) local gu,lu:
if k=0 then
 RETURN(X):
fi:

lu:=DerEqP(P,X,t):

if k=1 then
 RETURN(lu):
fi:



gu:=DerEqPk(P,X,t,k-1):
gu:=normal(t*normal(diff(gu,t)+diff(gu,X)*lu)):
gu:

end:

#AveAndMomsP(P,K): Inputs a losing (or symbolic) probability distribution P, finds a list of length K whose first entry is the expected duration until out of the game, the second is the standard deviation, followed by the 3rd through the
#K-th scaled moments. Try:
#AveAndMomsP([[1,4/10],[0,1/10],[-1,1/2]],4);
AveAndMomsP:=proc(P,K) local lu,X,t,mu,m2,sd,T,k,ku,i,T1,j:
lu:=DerEqP(P,X,t):

if subs({X=1,t=1},denom(lu))=0 then
# print(`Not yet implemented `):
 RETURN(FAIL):
fi:

mu:=subs({X=1,t=1},lu):

m2:=DerEqPk(P,X,t,2):


if subs({X=1,t=1},denom(m2))=0 then
 # print(`Not yet implemented `):
 RETURN(FAIL):
fi:

m2:=subs({X=1,t=1},m2):
sd:=sqrt(m2-mu^2):

T[1]:=mu:
T[2]:=m2:

for k from 3 to K do
ku:=DerEqPk(P,X,t,k):


if subs({X=1,t=1},denom(ku))=0 then
 #print(`Not yet implemented `):
 RETURN(FAIL):
fi:

ku:=subs({X=1,t=1},ku):
T[k]:=ku:
od:

T[0]:=1:
[seq(T[i],i=1..K)],sd:

for i from 3 to K do
 T1[i]:=normal(add(binomial(i,j)*T[j]*(-mu)^(i-j),j=0..i)):
 T1[i]:=T1[i]/sd^i:
od:

[mu,normal(sd),seq(normal(T1[i]),i=3..K)]:

end:




#AveAndMomsPc(st,P,K,L): Inputs a losing NUMERIC LOSING probability distribution P, a positive integer K, and a a (usually) large positive integer L
#STARTING with capital st
#outputs the a list of length K whose first entry is the CONDITIONAL expected duration until out of the game CONDITIONED that you exited after L steps
#the second is the CONDITIONAL standard deviation, followed by the 3rd through the
#K-th scaled moments. It also outputs the probability that you will exit in <=L rounds. Try:
#AveAndMomsPc(0,[[1,4/10],[0,1/10],[-1,1/2]],4,1000);
AveAndMomsPc:=proc(st,P,K,L) local f,s,t:
f:=GFpr(st,P,t,L):
s:=subs(t=1,f):
f:=f/s:
evalf(Stat(f,t,K)),evalf(s):

end:



#Paper1pS(k): An article about the algebraic equation satisfied by the probability generating function of life in a Casino where the
#probability of i (between -k and k) is p[i], and of course they all add-up to 1. try
#Paper1pS(1):

Paper1pS:=proc(k) local X,t,P,p,eq,n,t0,i,c0,eq1:
t0:=time():
P:=[seq([i,p[i]],i=-k..k)]:
eq:=EqP(P,X,t):
c0:=coeff(coeff(eq,X,1),t,0):
if c0=0 then
 RETURN(FAIL):
fi:
eq1:=expand(X-eq/c0):
eq:=collect(eq,X):
eq1:=collect(eq1,X):

print(``):
print(`The algebraic equation satisfied by the Probability Generating Function of Life in a Casion where the stakes are from`, -k, ` to `, k):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(``):
print(`Theorem:  Let a(n) be probability that you exit after EXACTLY n rounds in a casino where in each turn`):

for i from 1 to k do
print(`The probability that you win`, i, `dollars  is`, p[i]):
od:

for i from 1 to k do
print(`The probability that you lose`, i, `dollars  is`, p[-i]):
od:

print(`and the probability that you neither win nor lose is`, p[0]):

print(`Of course the probabilities have to add-up to 1, in other wrods`):
print(``):
print(add(p[i],i=-k..k)=1):
print(``):
print(`We also assume that the Casino stays in business, in other words, the expected gain in one round is negative:`):
print(``):
print(add(i*p[i],i=-k..k)<1):
print(``):

print(``):
print(`Let X(t) be the ordinary generating function of that sequence, in other words`):
print(``):
print(X(t)=Sum(a(n)*t^n,n=0..infinity)):
print(``):
print(`X(t) satisfies the algebraic equation `):
print(``):
print(subs(X=X(t),eq)=0):
print(``):
print(`and in Maple notation `):
print(``):
lprint(subs(X=X(t),eq)=0):
print(``):
print(`Or more usefully (for computing many terms)`):
print(``):
print(X(t)=subs(X=X(t),eq1)):
print(``):
print(`and in Maple notation `):
print(``):
lprint(X(t)=subs(X=X(t),eq1)):
print(``):

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

end:





#Paper1p(P,K1,K2,K3): An article about the algebraic equation satisfied by the probability generating function of life in a Casino where the
#you are given a finite probability distribution P. It also gives you either the exact (if possible), or conditioned upon exiting after L steps
#and does some confirming simulations with K3 tries. Try:
#Paper1p([[1,1/3],[0,1/6],[-1,1/2]],4,1000,10000);
Paper1p:=proc(P,K1,K2,K3) local X,t,eq,n,t0,i,c0,eq1,lu:
t0:=time():
if not (type(P,list) and type(K1,integer)  and type(K2,integer)  and type(K3,integer) ) then
 print(`bad input`):
 RETURN(FAIL):
fi:

if add(P[i][2],i=1..nops(P))<>1  then
 print(P, `is not a probability distribution `):
 RETURN(FAIL):
fi:


if min(seq(P[i][2],i=1..nops(P)))<0 then
 print(P, `is not a probability distribution `):
 RETURN(FAIL):
fi:



if add(P[i][1]*P[i][2],i=1..nops(P))>=0  then
 print(P, `is not a LOSING probability distribution, you will EXPECT to be in the Casino for ever, getting richer and richer `):
 RETURN(FAIL):
fi:

eq:=EqP(P,X,t):
c0:=coeff(coeff(eq,X,1),t,0):
if c0=0 then
 RETURN(FAIL):
fi:
eq1:=expand(X-eq/c0):
eq:=collect(eq,X):
eq1:=collect(eq1,X):

print(``):
print(`The algebraic equation satisfied by the Probability Generating Function of Life in a Casion where the stakes given by the Probability distribution`, P):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(``):
print(`Theorem:  Let a(n) be probability that you exit after EXACTLY n rounds in a casino where in each turn`):

for i from 1 to nops(P) do
if P[i][1]>0 then
print(`The probability that you win`, P[i][1], `dollars  is`, P[i][2]):
elif P[i][1]<0 then
print(`The probability that you lose`, -P[i][1], `dollars  is`, P[i][2]):
else
print(`The probability that you neither win nor lose is`, P[i][2]):
fi:
od:

print(``):
print(`Let X(t) be the ordinary generating function of that sequence, in other words`):
print(``):
print(X(t)=Sum(a(n)*t^n,n=0..infinity)):
print(``):
print(`X(t) satisfies the algebraic equation `):
print(``):
print(subs(X=X(t),eq)=0):
print(``):
print(`and in Maple notation `):
print(``):
lprint(subs(X=X(t),eq)=0):
print(``):
print(`Or more usefully (for computing many terms)`):
print(``):
print(X(t)=subs(X=X(t),eq1)):
print(``):
print(`and in Maple notation `):
print(``):
lprint(X(t)=subs(X=X(t),eq1)):
print(``):

lu:=AveAndMomsP(P,K1):

if lu<>FAIL then

 print(`The EXACT expected number of rounds  in the casino, until being in the red is`, lu[1], `or in decimals`, evalf(lu[1])):
print(``):
 print(`The EXACT standard-deviation of life in the casino is`, evalf(lu[2])):
for i from 3 to nops(lu) do
print(`The scaled`, i, `moment is`, evalf(lu[i])):
od:

else
lu:=AveAndMomsPc(0,P,K1,K2):
print(`The probability that you got broke  before`, K2, `rounds is`, lu[2]):
print(``):
lu:=evalf(lu[1]):
print(`assuming that this is the case `):

 print(`The Conditional  expected number of rounds  in the casino, until being in the red is`, lu[1]):
print(``):
 print(`The Conditonal standard-deviation of life in the casino is`, evalf(lu[2])):
for i from 3 to nops(lu) do
print(`The condintioanl scaled`, i, `moment is`, evalf(lu[i])):
od:


fi:


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

end:



#Paper2pS(k,K): An article about the algebraic equation satisfied by the probability generating function of life in a Casino where the
#probability of i (between  -1 and k) is p[i], and of course they all add-up to 1, and the expected gain is negative, it also gives you 
#explicit expressions for the expectation, standard-deviation, and the 3rd through the k-th moment in termss of the p[i]. try
#Paper2pS(1,4):

Paper2pS:=proc(k,K) local X,t,P,p,eq,n,t0,i,c0,eq1,lu:
t0:=time():
P:=[seq([i,p[i]],i=-1..k)]:
eq:=EqP(P,X,t):
c0:=coeff(coeff(eq,X,1),t,0):
if c0=0 then
 RETURN(FAIL):
fi:
eq1:=expand(X-eq/c0):
eq:=collect(eq,X):
eq1:=collect(eq1,X):

print(``):
print(`The algebraic equation satisfied by the Probability Generating Function of Life in a Casion where the stakes are from`, -1, ` to `, k):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(``):
print(`Theorem:  Let a(n) be probability that you exit after EXACTLY n rounds in a casino where in each turn`):

for i from 1 to k do
print(`The probability that you win`, i, `dollars  is`, p[i]):
od:

for i from -1 to -1 do
print(`The probability that you lose`, i, `dollars  is`, p[i]):
od:

print(`and the probability that you neither win nor lose is`, p[0]):

print(`Of course the probabilities have to add-up to 1, in other wrods`):
print(``):
print(add(p[i],i=-1..k)=1):
print(``):
print(`We also assume that the Casino stays in business, in other words, the expected gain in one round is negative:`):
print(``):
print(add(i*p[i],i=-1..k)<1):
print(``):

print(``):
print(`Let X(t) be the ordinary generating function of that sequence, in other words`):
print(``):
print(X(t)=Sum(a(n)*t^n,n=0..infinity)):
print(``):
print(`X(t) satisfies the algebraic equation `):
print(``):
print(subs(X=X(t),eq)=0):
print(``):
print(`and in Maple notation `):
print(``):
lprint(subs(X=X(t),eq)=0):
print(``):
print(`Or more usefully (for computing many terms)`):
print(``):
print(X(t)=subs(X=X(t),eq1)):
print(``):
print(`and in Maple notation `):
print(``):
lprint(X(t)=subs(X=X(t),eq1)):
print(``):

lu:=AveAndMomsP(P,K):

print(`The explicit expression for the  expected duration in the casino is`):
print(``):
print(lu[1]):
print(``):
print(`and in Maple notation `):
print(``):
print(``):
lprint(lu[1]):
print(``):


print(`The explicit expression for the  standard-deviation is`):
print(``):
print(lu[2]):
print(``):
print(`and in Maple notation `):
print(``):
print(``):
lprint(lu[2]):
print(``):


for i from 3 to K do
print(``):

print(`The explicit expression for the`, i, `-th scaled moment is`):
print(``):
print(lu[i]):
print(``):
print(`and in Maple notation `):
print(``):
print(``):
lprint(lu[i]):
print(``):

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

end:



#AveSD(P): The exact average and standard-deviation  of duration with NUMERIC prob. dist. P.
#Try:
#AveSD([[1,1/4],[0,1/4],[-1,1/2]]);
AveSD:=proc(P) local eq,i,mu,f,m2,X,t,lu:


if {seq(type(P[i][2],numeric),i=1..nops(P))}<>{true} then
 print(P, `had to be numeric`):
 RETURN(FAIL):
fi:


if add(P[i][2],i=1..nops(P))<>1  then
 print(P, `is not a probability distribution `):
 RETURN(FAIL):
fi:


if min(seq(P[i][2],i=1..nops(P)))<0 then
 print(P, `is not a probability distribution `):
 RETURN(FAIL):
fi:



if add(P[i][1]*P[i][2],i=1..nops(P))>=0  then
 print(P, `is not a LOSING probability distribution, you will EXPECT to be in the Casino for ever, getting richer and richer `):
 RETURN(FAIL):
fi:



mu:=AveAndMomsP(P,2):
if mu<>FAIL then
 RETURN(mu):
fi:

eq:=EqP(P,X,t):
eq:=expand(subs({X=X+1,t=t+1},eq)):
f:=solve(eq,X);

f:=taylor(f,t=0,5):
mu:=coeff(f,t,1):
m2:=coeff(f,t,2)*2+mu:

lu:=AveAndMomsPc(0,P,2,400)[1]:
if evalf(abs(lu[1]-mu))>1/10^4 then
 print(`something is fishy`):
 RETURN(FAIL):
fi:

[mu,sqrt(m2-mu^2)]:

end:






#MakeSysT1(f,g,t,S,f[0,0],X): converts MakeSysT(f,g,t,S) to a transoformation format, [LIST,LIST]
MakeSysT1:=proc(f,g,t,S,st,X) local gu,eq,var,var1,eq1,i,T,Eqs:

gu:=MakeSysT(f,g,t,S):

eq:=gu[1]:
var:=gu[2]:

if not member(st,var) then
 RETURN(FAIL):
fi:

var1:= var minus {st}:
var:=[st,op(convert(var1,list))]:


for eq1 in eq do
 T[op(1,eq1)]:=op(2,eq1):
od:

Eqs:=[seq(T[var[i]],i=1..nops(var))]:

for i from 1 to nops(var) do
T[var[i]]:=X[i]:
od:

Eqs:=subs({seq(var[i]=T[var[i]],i=1..nops(var))},Eqs):
end:







#Trunc1(f,t,K): truncates the polynomial f in t up to degree K
Trunc1:=proc(f,t,K) local i:
add(coeff(f,t,i)*t^i,i=0..K):
end:


#SeqFromSysEq(Eqs,t,X,K): Given a system of  equations in terms of X[1],X[2],... and a variable T
#finds the first K terms of X[1]. Try:
#gu:=MakeSysT1(f,g,t,{1,-1,0},f[0,0],X): SeqFromSysEq(gu,t,X,10);
SeqFromSysEq:=proc(Eqs,t,X,K) local var1,var2,var3,k,i:
k:=nops(Eqs):
var1:=[1,0$(k-1)]:

var2:=expand(subs({seq(X[i]=var1[i],i=1..k)},Eqs)):
var2:=[seq(Trunc1(var2[i],t,K),i=1..nops(var2))]:

while var1<>var2 do
var3:=expand(subs({seq(X[i]=var2[i],i=1..k)},Eqs)):
var3:=[seq(Trunc1(var3[i],t,K),i=1..nops(var2))]:
var1:=var2:
var2:=var3:
od:

var2:

end:



#SeqSee(S,K): Same as SeqS(S,K) and SeqSe(S,K), but using the System of equation satisfied by the generating function obtained via EqGFt(S,X,t). Try:
#SeqSee({1,0,-1},30);
SeqSee:=proc(S,K) local gu,f,g,X,i,t:

gu:=MakeSysT1(f,g,t,S,f[0,0],X): 

gu:=SeqFromSysEq(gu,t,X,K)[1];

[seq(coeff(gu,t,i),i=0..K)]:


end:

###FROM qGDW.txt



#CheckEq(eq,X,t,L): Given an equation eq in terms of X and t and a sequence L starting at n=0, checks whether
#add(L[i]*t^i,i=1..nops(L)) satisfies it. Try
#CheckEq(X-1-t*X^2,X,t,[seq(binomial(2*i,i)/(i+1),i=0..50)]);
CheckEq:=proc(eq,X,t,L) local X1,i:
X1:=add(L[i]*t^(i-1),i=1..nops(L)):
if ldegree(expand(subs(X=X1,eq)),t)>nops(L)-2 then
 true:
else
 false:
fi:
end:

qnmwd:=proc(S,n,m,q) local W,S2,s,W1: 
	option remember:
	W := 0:
	S2 := {seq(s[2],s in S)}:
	if n=0 then 
		RETURN(1):
	elif n=1 then 
		if member(m,S2)=true and 0<=m then 
			RETURN(q^m) else RETURN(0): 
		fi:
	else 
		for s in S2 do 
			if 0<=m - s then :
				W1 := qnmwd(S,n - 1,m - s,q):  
				if W1<>0 then 
					W := expand(W+q^m*W1):
				fi:
			fi:
		od:
	fi:	
W:
end:


qnwd:=proc(S,n,q) option remember:
	qnmwd(S,n,0,q):
end:

qnwdK:=proc(S,K,q) local n: 
	[seq(qnwd(S,n,q),n=1..K)]:
end:


#qSeqS(S,q,K): The first K terms, starting with n=0 of the weight-enumerator, according to q^area,  of generalized Dyck path with set steps S .
#Try:
#qSeqS({1,0,-1},q,10);
qSeqS:=proc(S,q,K) local s:
[1,op(qnwdK({seq([1,s], s in S)},K,q))] :
end:


#qMakeEqFt(f, g, t, q,a, b, S): Sets up a functional equation for the quantity f[a,b](t,q), the set of walks using the set of steps {[1,s], s in S}, that start at [0,a] and end where y=b
#and that are WEAKLY above the x-axis. Try:
#qMakeEqFt(f,g,t,q,0,0,{1,0,-1});
qMakeEqFt:= proc(f, g, t, q,a, b, S):
	if 0 < a and 0 < b then 
		RETURN(f[a, b](t) -( f[0, b](t)*g[a, 0](t) + f[a - 1, b - 1](q*t)), {f[a,b](t), f[0, b](t), f[a - 1, b - 1](q*t), g[a, 0](t)}, {f[a,b],f[0,b],f[a-1,b-1]}   ):


	elif a = 0 and b > 0 then 
		RETURN(f[a, b](t) -( f[0, 0](t)*g[0, b](t)), {f[a,b](t),f[0, 0](t), g[0, b](t)}, {f[a,b],f[0,0],g[0,b]} ): 

	elif 0 < a and b = 0 then 
		RETURN(f[a, b](t) -( g[a, 0](t)*f[0, 0](t)), {f[a,b](t),f[0, 0](t), g[a, 0](t)}, {f[a,b],f[0,0], g[a,0]}):

 	elif a = 0 and b = 0 then 
		if not member(0, S) then 
			RETURN(f[0, 0](t) -( f[0, 0](t)*g[0, 0](t) + 1), {f[0, 0](t), g[0, 0](t)},{f[0,0],g[0,0]}):

		else 
			RETURN(f[0, 0](t) -(f[0, 0](t)*g[0, 0](t) + f[0, 0](t)*t + 1), {f[0, 0](t), g[0, 0](t)},{f[0,0],g[0,0]}):
fi:
fi:
end:



#qMakeEqGt(f, g, t, q,a, b, S): Sets up a functional equation for the quantity g[a,b](t,q), the set of walks using the set of steps {[1,s], s in S}, that start at [0,a] and end where y=b
#and that are STRICTLY above the x-axis (expect possibly at the endpoints). Try:
#qMakeEqGt(f,g,t,q,0,0,{1,0,-1});

qMakeEqGt := proc(f, g, t,q, a, b, S) local gu, s, P, N, s1, s2, lu,mu:
	P := {}:
	 N := {}:

	 for s in S do 
		if 0 < s then 
			P := P union {s}:
		elif s < 0 then 
			N := N union {s}:
		fi:
	od:

	if 0 < a and 0 < b then 
		print(`Something is wrong`): 
		RETURN(FAIL): 
	fi:

	gu := 0: lu := {g[a,b](t)}: mu:={g[a,b]}:

	 if a = 0 and 0 < b then 
		for s in P do 
			lu := lu union {f[a + s - 1, b - 1](t), f[a + s - 1, b - 1](q*t) }:
                     mu := mu union {f[a + s - 1, b - 1] }:
			gu := gu + t*q^(1/2*s)*f[a + s - 1, b - 1](q*t):
		od:
 
		 RETURN(g[a, b](t)-gu, lu,mu):
	fi:

	if 0 < a and b = 0 then 
		for s in N do 
			lu := lu union {f[a - 1, b - s - 1](t), f[a - 1, b - s - 1](q*t) }:
			mu := mu union {f[a - 1, b - s - 1] }:
			gu := gu + t*q^(-1/2*s) * f[a - 1, b - s - 1](q*t):
		od:

		RETURN(g[a, b] (t)- gu, lu,mu):
	fi:

	if a = 0 and b = 0 then 
		for s1 in P do 
			for s2 in N do 
				lu := lu union {f[a + s1 - 1, b - s2 - 1](t), f[a + s1 - 1, b - s2 - 1](q*t) }:
                                mu := mu union {f[a + s1 - 1, b - s2 - 1] }:
				gu := gu + t*q^(a + 1/2*s1)*t*q^(-1/2*s2)*f[a + s1 - 1, b - s2 - 1](q*t):
			od:
		od:
 		
		RETURN(g[0, 0](t) - gu, lu,mu):
	fi:
end:

#qMakeSyst(f,g,t,q,S): inputs symbols f, g, for functions, and t,q, for variables, and a set of (positive and negative integers) set up a system of
#functional equations for the quantities f[a,b](t,q) and g[a,b](t,q), where 
#f[a,b](t,q) (respectively g[a,b](q)) is the weight-enumerator according to the weight
#t^length(w)*q^area(w) of all paths using the steps {[1,s], s in S} starting at [0,a] and ending at [n,b] and WEAKLY (resp. STRICTLY, except possibly at the endpoints) above the x-axis.
#For example, the system for Motzkin paths is gotten by typing:
#qMakeSys(f,g,t,q,{1,0,-1});
qMakeSyst:= proc(f, g, t, q, S) local eq, StillToDo, AlreadyDone, guy, ku,VARS:
	eq := {}:
	StillToDo := {f[0, 0]}:
	AlreadyDone := {}:
       VARS:={}:
	while StillToDo <> {} do 
		guy := StillToDo[1]:

		if op(0, guy) = f then 
			ku := qMakeEqFt(f, g, t, q, op(1, guy), op(2, guy), S):
			AlreadyDone := AlreadyDone union {guy}:
			StillToDo := (StillToDo union ku[3]) minus AlreadyDone:
                        VARS:=VARS union ku[2]:
			eq := eq union {ku[1]}:
		elif  op(0, guy) = g then 
			ku := qMakeEqGt(f, g, t, q, op(1, guy), op(2, guy), S):
			AlreadyDone := AlreadyDone union {guy}:
 			StillToDo := (StillToDo union ku[3]) minus AlreadyDone:
                        VARS:=VARS union ku[2]:
		 	eq := eq union {ku[1]}:
 		else 
			print(`Something is wrong`):
			 RETURN(FAIL): 
		fi:
	od: 

	eq, VARS:
end:


#qEqGFt(S,X,t): Inputs a set of positive and negative integers S, outputs an algebraic equation for X=X(t) the
#generating function for the sun-of-the-areas under all generalized Dyck paths with steps {[1,s], s in S}. For the
#classical Dyck case do:
#qEqGFt({1,-1},X,t);

qEqGFt:=proc(S,X,t) local q,f,g,sys, eq, var,eq1,v,EQ0,EQ1,EQ2,EQ,VAR,i,Z,lu,L,lu1,aj,S1,s:
aj:=igcd(op(S)):

if aj>1 then
 S1:={seq(s/aj,s in S)}:
 RETURN(numer(subs(X=X/aj,qEqGFt(S1,X,t)))):
fi:
	sys:=qMakeSyst(f,g,t,q,S):
	eq:=sys[1]:
	var:=sys[2]: 

	eq1:=expand(subs( {seq(v = op(0, op(0, v))[op(1 .. 2, op(0, v)), 0](op(1, v)) + (q - 1)*op(0, op(0, v))[op(1 .. 2, op(0, v)), 1](op(1, v)), v in var)},eq) 
	): 
	var:=subs(q=1,var):  
	VAR:={seq(op(0,op(0,v))[op(1..2,op(0,v)),0](op(1,v)),v in var),seq(op(0,op(0,v))[op(1..2,op(0,v)),1](op(1,v)),v in var),seq(diff(op(0,op(0,v))[op(1..2,op(0,v)),0](op(1,v)),op(1,v)),v in var)}:  

	EQ0:=subs(q=1,eq1):
	EQ1:=subs(q=1,diff(eq1,q)):
	EQ2:=subs(q=1,diff(eq1,t)):
	EQ:=EQ0 union EQ1 union EQ2: 
	EQ:=subs({seq(D(op(0,op(0,v))[op(1..2,op(0,v)),0]) (op(1,v))=diff(op(0,op(0,v))[op(1..2,op(0,v)),0](op(1,v)),op(1,v)),v in var)},EQ):
	VAR:=[op(VAR minus {f[0,0,1](t)}),f[0,0,1](t)]:
	EQ:=subs({seq(VAR[i]=Z[i],i=1..nops(VAR))},EQ): 
        #print(EQ):
	VAR:=subs({seq(VAR[i]=Z[i],i=1..nops(VAR))},VAR): 
	lu:=subs(Z[nops(VAR)] = X, Groebner[Basis](EQ, plex(op(VAR)))[1]):
	lu:=add(factor(coeff(lu,X,i))*X^i,i=0..degree(lu,X)):
     lu:

    L:=qSeqS(S,q,50):

    L:=[seq(subs(q=1,diff(L[i],q)),i=1..nops(L))]:

 if not CheckEq(lu,X,t,L) then
  lprint(lu):
 print(`it did not work out`):
 RETURN(FAIL):
 else

lu:=factor(lu):

if not type(lu,`*`) then
lu:=add(factor(coeff(lu,X,i))*X^i,i=0..degree(lu,X)):
 RETURN(lu):
fi:


for i from 1 to nops(lu) do
lu1:=op(i,lu):
if CheckEq(lu1,X,t,L) then
lu1:=add(factor(coeff(lu1,X,i))*X^i,i=0..degree(lu1,X)):
 RETURN(lu1):
fi:
od:
eq:

lu:=add(factor(coeff(lu,X,i))*X^i,i=0..degree(lu,X)):
 RETURN(lu):
 fi:
end:




#Paper1A(S): An article about the enumerating sequence of words in the alphabet S that add up to 0 and whose partial sums are never neagtive, and about the sum of the areas/ Try
#Paper1A({1,0,-1}):

Paper1A:=proc(S) local X,t,Y,t0,a,b,i,n,q,s,L,eq:
t0:=time():
print(``):
print(`Enumerating generalized Dyck words in the alphabet`, S, ` as well as the Sum of the Areas Under Them`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
eq:=EqGFt(S,X,t):
eq:=add(factor(coeff(eq,Y,i))*X^i,i=0..degree(eq,Y)):
print(``):
print(`Theorem:  Let a(n) be the number of words of length n in the alphabet`, S, ` that sum-up to 0 and whose partial sums are never negative, in other words generalized Dyck words with alphabet`, S, `rather than {1,-1} `):
print(``):
print(`Let X(t) be the ordinary generating function of that sequence, in other words`):
print(``):
print(X(t)=Sum(a(n)*t^n,n=0..infinity)):
print(``):
print(`X(t) satisfies the algebraic equation `):
print(``):
print(subs(X=X(t),eq)=0):
print(``):
print(`and in Maple notation `):
print(``):
lprint(subs(X=X(t),eq)=0):
print(``):
print(`For the sake of the OEIS here are the first 41 terms, starting with n=0`):
print(``):
lprint(SeqS(S,40)):
print(``):
print(`Using the algebraic equation, we can get many more terms. Here are the first 301 terms, starting at n=0`):
print(``):
lprint(AlgToSeq(eq,X,t,300)):
print(``):


eq:=qEqGFt(S,Y,t):
eq:=add(factor(coeff(eq,Y,i))*Y^i,i=0..degree(eq,Y)):
print(``):
print(`Theorem:  Let b(n) be the sum of the areas under generalized Dyck paths of length n in the set of steps`, {seq([1,s], s in S)}):
print(``):
print(`Let Y(t) be the ordinary generating function of that sequence, in other words`):
print(``):
print(Y(t)=Sum(b(n)*t^n,n=0..infinity)):
print(``):
print(`Y(t) satisfies the algebraic equation `):
print(``):
print(subs(Y=Y(t),eq)=0):
print(``):
print(`and in Maple notation `):
print(``):
lprint(subs(Y=Y(t),eq)=0):
print(``):
print(`For the sake of the OEIS here are the first 41 terms, starting with n=0`):
print(``):
L:=qSeqS(S,q,40):
lprint([seq(subs(q=1,diff(L[i],q)),i=1..nops(L))]):
print(``):

print(``):
print(`Using the algebraic equation, we can get many more terms. Here are the first 301 terms, starting at n=0`):
print(``):
lprint(AlgToSeq(eq,Y,t,300)):
print(``):

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



#AlgToSeq(eq,X,t,K): Given an algebraic equation in terms of X=X(t) and t, finds the first K terms starting at n=0. Try:
#AlgToSeq(X-1-t*X^2,X,t,30);
AlgToSeq:=proc(eq,X,t,K) local gu,eq1,i,c,gu1,gu2:
c:=coeff(coeff(eq,t,0),X,1):

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

eq1:=eq-c*X:
eq1:=expand(-eq1/c):

gu:=1:

gu1:=expand(subs(X=gu,eq1)):
gu1:=add(coeff(gu1,t,i)*t^i,i=0..K):

while gu<>gu1 do

gu2:=expand(subs(X=gu1,eq1)):
gu2:=add(coeff(gu2,t,i)*t^i,i=0..K):
gu:=gu1:
gu1:=gu2:
od:

[seq(coeff(gu2,t,i),i=0..K)]:

end:




#Paper11A(S,c): An article about the enumerating sequence of words in the alphabet S that add up to 0 and whose partial sums are never neagtive, and about the sum of the areas with a numbered theorem c. Try
#Paper11A({1,0,-1},3):

Paper11A:=proc(S,c) local X,t,Y,t0,a,b,i,n,s,eq,eq1,L,L1,R:
t0:=time():
print(``):
eq:=EqGFt(S,X,t):
eq1:=qEqGFt(S,Y,t):
eq:=add(factor(coeff(eq,X,i))*X^i,i=0..degree(eq,X)):
eq1:=add(factor(coeff(eq1,Y,i))*Y^i,i=0..degree(eq1,Y)):

print(``):
print(`Theorem Number`, c,  `Let a(n) be the number of generalized  Dyck paths of length n in the set of steps`, {seq([1,s], s in S)}):
print(`also Let b(n) be the sum of the areas under these generalized Dyck paths of length n in the set of steps`):
print(``):
print(`Let X(t), Y(t)  be the ordinary generating functions of the sequences a(n), b(n), in other words`):
print(``):
print(X(t)=Sum(a(n)*t^n,n=0..infinity)):
print(``):
print(Y(t)=Sum(b(n)*t^n,n=0..infinity)):
print(``):
print(`X(t),satisfies the algebraic equation `):
print(``):
print(subs(X=X(t),eq)=0):
print(``):
print(`and in Maple notation `):
print(``):
lprint(subs(X=X(t),eq)=0):
print(``):

print(`Y(t),satisfies the algebraic equation `):
print(``):
print(subs(Y=Y(t),eq1)=0):
print(``):
print(`and in Maple notation `):
print(``):
lprint(subs(Y=Y(t),eq1)=0):
print(``):

print(``):


L:=AlgToSeq(eq,X,t,200):
L1:=AlgToSeq(eq1,Y,t,200):
print(`Using these algebraic equations, we can get many more terms. Here are the first 101 terms of the  enuerating sequence, a(n), starting at n=0, are`):
print(``):
lprint([op(1..101,L)]):
print(``):

print(`The first 101 terms of the  sum of areas sequence, b(n), starting at n=0, are`):
print(``):
lprint([op(1..101,L1)]):
print(``):

print(`Finally, here are the average areas for n from 190 to 200 divided by n^(3/2)`):

R:=[]:
for i from 190 to 200 do
 if L[i+1]<>0 then
R:=[op(R),[i,evalf(L1[i+1]/L[i+1]/i^(3/2))]]:
fi:
od:


print(``):
lprint(R):

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




#BookA(M): Inputs a positive integer M and outputs a book about enumearing generalized Dyck path with all possible alphabets consisting of both negative and positive integes and 0. It also
#BookA(1);
BookA:=proc(M) local NEG,POS,c,i,gu1,gu2,gu11,gu21,t0,S:
t0:=time():
print(``):
print(`Enumerating Generalized Dyck paths and The Sum of Areas with alphabets consisting of integers from`, -M, ` to `, M):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
c:=0:
NEG:={seq(-i,i=1..M)}:
POS:={seq(i,i=1..M)}:
gu1:=powerset(NEG) minus {{}}:
gu2:=powerset(POS) minus {{}}:
print(``):
print(`-------------------------------------------------------------`):
print(``):

for gu11 in gu1 do
for gu21 in gu2 do
 S:=gu11 union gu21:

if igcd(op(S))=1 then

c:=c+1:

Paper11A(S,c):
print(``):
print(`-------------------------------------------------------------`):
print(``):
S:=S union {0}:
c:=c+1:
Paper11A(S,c):
fi:
od:
od:

print(``):
print(`--------------------------------`):
print(``):
print(`This ends this book that took`, time()-t0, `seconds to create `):
print(``):
end:






#Paper2A(S,MaxC): An article about the sequence giving the number, and sum of areas, of generalized Dyck path with set of steps S and tries to find recurrences
#for them if their complexity is <=MaxC. Otherwise it just gives the beginning of these sequences. Try:
#Paper2A({1,0,-1},10):

Paper2A:=proc(S,MaxC) local t0,K,L,L0,L1,L2,D1,ope0,ope1,ope2,i,q,n,N,a,s,ope0N,ope1N,ope2N,gu:
t0:=time():
print(``):
print(`Searching for Recurrences Satisfied by the Sequences Enumerating  generalized Dyck words in the alphabet`, S, ` as well as the Sum of the Areas Under Them Sequence of complexity <=`, MaxC):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

K:=max(seq((2+D1)*(1+MaxC-D1),D1=0..MaxC))+5:

L:=qSeqS(S,q,K):
L:=[op(2..nops(L),L)]:
L0:=[seq(subs(q=1,L[i]),i=1..nops(L))]:
L1:=[seq(subs(q=1,diff(L[i],q)),i=1..nops(L))]:
L2:=[seq(subs(q=1, diff(q*diff(L[i],q),q)),i=1..nops(L))]:



ope0:=Findrec(L0,n,N,MaxC):

print(`Let a(n) be the NUMBER of generalized Dyck path with set of steps`, S):
print(`(i.e. walks from [0,0] to [n,0] with atomic steps`, {seq([1,s], s in S)}, `that NEVER dip below the x-axis`):


if ope0=FAIL then
print(`There is no homog. linear recurrence with ORDER+DEGREE<=`, MaxC, ` satisfied by a(n), so here are the first `, K, `terms. `):
print(``):
print(L0):
else
ope0N:=expand(1-expand(subs(n=n-degree(ope0,N),ope0)/N^degree(ope0,N))):


ope0N:=add(factor(coeff(ope0N,N,-i))*N^(-i),i=1..-ldegree(ope0N,N)):

print(`The sequence a(n) satisfies the following linear recurrence equation`):

print(``):
print(a(n)=add(coeff(ope0N,N,-i)*a(n-i),i=1..-ldegree(ope0N,N))):
print(``):
print(`and in Maple notation`):
print(``):
lprint(a(n)=add(coeff(ope0N,N,-i)*a(n-i),i=1..-ldegree(ope0N,N))):
print(``):

print(`Subject to the initial conditions`):
print(``):
lprint(seq(a(i)=L0[i],i=1..degree(ope0,N))):
print(``):
L0:=SeqFromRec(ope0,n,N,[op(1..degree(ope0,N),L0)],1010):
L0:=[op(1001..1010,L0)]:
print(``):
print(`Here are the terms of a(n) from n=1001 to n=1010 `):
print(``):
lprint(L0):
fi:



print(`----------------------------------------------`):
print(``):
ope1:=Findrec(L1,n,N,MaxC):

print(`Let b(n) be the SUM OF THE AREAS UNDER  generalized Dyck path with set of steps`, S):
print(`(i.e. walks from [0,0] to [n,0] with atomic steps`, {seq([1,s], s in S)}, `that NEVER dip below the x-axis`):


if ope1=FAIL then
print(`There is no homog. linear recurrence with ORDER+DEGREE<=`, MaxC, ` satisfied by b(n), so here are the first `, K, `terms. `):
print(``):
print(L1):
else
ope1N:=expand(1-expand(subs(n=n-degree(ope1,N),ope1)/N^degree(ope1,N))):


ope1N:=add(factor(coeff(ope1N,N,-i))*N^(-i),i=1..-ldegree(ope1N,N)):

print(`The sequence b(n) satisfies the following linear recurrence equation`):

print(``):
print(b(n)=add(coeff(ope1N,N,-i)*b(n-i),i=1..-ldegree(ope1N,N))):
print(``):
print(`and in Maple notation`):
print(``):
lprint(b(n)=add(coeff(ope1N,N,-i)*b(n-i),i=1..-ldegree(ope1N,N))):
print(``):

print(`Subject to the initial conditions`):
print(``):
lprint(seq(b(i)=L1[i],i=1..degree(ope0,N))):
print(``):
L1:=SeqFromRec(ope1,n,N,[op(1..degree(ope1,N),L1)],1010):
L1:=[op(1001..1010,L1)]:
print(``):
print(`Here are the terms of b(n) from n=1001 to n=1010 `):
print(``):
lprint(L1):
fi:




print(`----------------------------------------------`):
print(``):
ope2:=Findrec(L2,n,N,MaxC):

print(`Let c(n) be the SUM OF THE SQUARED-AREAS UNDER  generalized Dyck path with set of steps`, S):
print(`(i.e. walks from [0,0] to [n,0] with atomic steps`, {seq([1,s], s in S)}, `that NEVER dip below the x-axis`):


if ope2=FAIL then
print(`There is no homog. linear recurrence with ORDER+DEGREE<=`, MaxC, ` satisfied by c(n), so here are the first `, K, `terms. `):
print(``):
print(L2):
else
ope2N:=expand(1-expand(subs(n=n-degree(ope2,N),ope2)/N^degree(ope2,N))):


ope2N:=add(factor(coeff(ope2N,N,-i))*N^(-i),i=1..-ldegree(ope2N,N)):

print(`The sequence c(n) satisfies the following linear recurrence equation`):

print(``):
print(c(n)=add(coeff(ope2N,N,-i)*c(n-i),i=1..-ldegree(ope2N,N))):
print(``):
print(`and in Maple notation`):
print(``):
lprint(c(n)=add(coeff(ope2N,N,-i)*c(n-i),i=1..-ldegree(ope2N,N))):
print(``):

print(`Subject to the initial conditions`):
print(``):
lprint(seq(c(i)=L2[i],i=1..degree(ope2,N))):
print(``):
L2:=SeqFromRec(ope2,n,N,[op(1..degree(ope2,N),L2)],1010):
L2:=[op(1001..1010,L2)]:
print(``):
print(`Here are the terms of c(n) from n=1001 to n=1010 `):
print(``):
lprint(L2):
fi:

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

if ope0<>FAIL and ope1<>FAIL then
print(`The average areas divided by n^(3/2) for n from 1000 to 1010 are`):
gu:=evalf([seq(L1[i]/L0[i]/(1000+i)^(3/2),i=1..10)]):
print(``):
print(gu):
print(``):
fi:




if ope0<>FAIL and ope1<>FAIL and ope2<>FAIL then

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

print(`The variance of the area divided by n^3 for n from 1000 to 1010 are`):
gu:=evalf([seq( (L2[i]/L0[i]-(L1[i]/L0[i])^2) /(1000+i)^3,i=1..10)]):
print(``):
print(gu):
print(``):
fi:



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






#Paper2Ac(S,MaxC,c): An article about the sequence giving the number, and sum of areas, of generalized Dyck path with set of steps S and tries to find recurrences
#for them if their complexity is <=MaxC. Otherwise it just gives the beginning of these sequences. Try:
#Paper2Ac({1,0,-1},10,1):

Paper2Ac:=proc(S,MaxC,c1) local t0,K,L,L0,L1,L2,D1,ope0,ope1,ope2,i,q,n,N,a,s,ope0N,ope1N,ope2N,gu:
t0:=time():
print(``):
print(`Section Number`, c1):
print(``):
print(`Regarding the set of steps`, S):

K:=max(seq((2+D1)*(1+MaxC-D1),D1=0..MaxC))+5:

L:=qSeqS(S,q,K):
L:=[op(2..nops(L),L)]:
L0:=[seq(subs(q=1,L[i]),i=1..nops(L))]:
L1:=[seq(subs(q=1,diff(L[i],q)),i=1..nops(L))]:
L2:=[seq(subs(q=1, diff(q*diff(L[i],q),q)),i=1..nops(L))]:



ope0:=Findrec(L0,n,N,MaxC):

print(`Let a(n) be the NUMBER of generalized Dyck path with set of steps`, S):
print(`(i.e. walks from [0,0] to [n,0] with atomic steps`, {seq([1,s], s in S)}, `that NEVER dip below the x-axis`):


if ope0=FAIL then
print(`There is no homog. linear recurrence with ORDER+DEGREE<=`, MaxC, ` satisfied by a(n), so here are the first `, K, `terms. `):
print(``):
print(L0):
else
ope0N:=expand(1-expand(subs(n=n-degree(ope0,N),ope0)/N^degree(ope0,N))):


ope0N:=add(factor(coeff(ope0N,N,-i))*N^(-i),i=1..-ldegree(ope0N,N)):

print(`The sequence a(n) satisfies the following linear recurrence equation`):

print(``):
print(a(n)=add(coeff(ope0N,N,-i)*a(n-i),i=1..-ldegree(ope0N,N))):
print(``):
print(`and in Maple notation`):
print(``):
lprint(a(n)=add(coeff(ope0N,N,-i)*a(n-i),i=1..-ldegree(ope0N,N))):
print(``):

print(`Subject to the initial conditions`):
print(``):
lprint(seq(a(i)=L0[i],i=1..degree(ope0,N))):
print(``):
L0:=SeqFromRec(ope0,n,N,[op(1..degree(ope0,N),L0)],1010):
L0:=[op(1001..1010,L0)]:
print(``):
print(`Here are the terms of a(n) from n=1001 to n=1010 `):
print(``):
lprint(L0):
fi:



print(`----------------------------------------------`):
print(``):
ope1:=Findrec(L1,n,N,MaxC):

print(`Let b(n) be the SUM OF THE AREAS UNDER  generalized Dyck path with set of steps`, S):
print(`(i.e. walks from [0,0] to [n,0] with atomic steps`, {seq([1,s], s in S)}, `that NEVER dip below the x-axis`):


if ope1=FAIL then
print(`There is no homog. linear recurrence with ORDER+DEGREE<=`, MaxC, ` satisfied by b(n), so here are the first `, K, `terms. `):
print(``):
print(L1):
else
ope1N:=expand(1-expand(subs(n=n-degree(ope1,N),ope1)/N^degree(ope1,N))):


ope1N:=add(factor(coeff(ope1N,N,-i))*N^(-i),i=1..-ldegree(ope1N,N)):

print(`The sequence b(n) satisfies the following linear recurrence equation`):

print(``):
print(b(n)=add(coeff(ope1N,N,-i)*b(n-i),i=1..-ldegree(ope1N,N))):
print(``):
print(`and in Maple notation`):
print(``):
lprint(b(n)=add(coeff(ope1N,N,-i)*b(n-i),i=1..-ldegree(ope1N,N))):
print(``):

print(`Subject to the initial conditions`):
print(``):
lprint(seq(b(i)=L1[i],i=1..degree(ope0,N))):
print(``):
L1:=SeqFromRec(ope1,n,N,[op(1..degree(ope1,N),L1)],1010):
L1:=[op(1001..1010,L1)]:
print(``):
print(`Here are the terms of b(n) from n=1001 to n=1010 `):
print(``):
lprint(L1):
fi:




print(`----------------------------------------------`):
print(``):
ope2:=Findrec(L2,n,N,MaxC):

print(`Let c(n) be the SUM OF THE SQUARED-AREAS UNDER  generalized Dyck path with set of steps`, S):
print(`(i.e. walks from [0,0] to [n,0] with atomic steps`, {seq([1,s], s in S)}, `that NEVER dip below the x-axis`):


if ope2=FAIL then
print(`There is no homog. linear recurrence with ORDER+DEGREE<=`, MaxC, ` satisfied by c(n), so here are the first `, K, `terms. `):
print(``):
print(L2):
else
ope2N:=expand(1-expand(subs(n=n-degree(ope2,N),ope2)/N^degree(ope2,N))):


ope2N:=add(factor(coeff(ope2N,N,-i))*N^(-i),i=1..-ldegree(ope2N,N)):

print(`The sequence c(n) satisfies the following linear recurrence equation`):

print(``):
print(c(n)=add(coeff(ope2N,N,-i)*c(n-i),i=1..-ldegree(ope2N,N))):
print(``):
print(`and in Maple notation`):
print(``):
lprint(c(n)=add(coeff(ope2N,N,-i)*c(n-i),i=1..-ldegree(ope2N,N))):
print(``):

print(`Subject to the initial conditions`):
print(``):
lprint(seq(c(i)=L2[i],i=1..degree(ope2,N))):
print(``):
L2:=SeqFromRec(ope2,n,N,[op(1..degree(ope2,N),L2)],1010):
L2:=[op(1001..1010,L2)]:
print(``):
print(`Here are the terms of c(n) from n=1001 to n=1010 `):
print(``):
lprint(L2):
fi:

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

if ope0<>FAIL and ope1<>FAIL then
print(`The average areas divided by n^(3/2) for n from 1000 to 1010 are`):
gu:=evalf([seq(L1[i]/L0[i]/(1000+i)^(3/2),i=1..10)]):
print(``):
print(gu):
print(``):
fi:




if ope0<>FAIL and ope1<>FAIL and ope2<>FAIL then

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

print(`The variance of the area divided by n^3 for n from 1000 to 1010 are`):
gu:=evalf([seq( (L2[i]/L0[i]-(L1[i]/L0[i])^2) /(1000+i)^3,i=1..10)]):
print(``):
print(gu):
print(``):
fi:



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



#BookA2(M,MaxC): Inputs a positive integer M and outputs a book about enumearing generalized Dyck path as well as the sum of the areas and area-squared
#with all possible alphabets consisting of both negative and positive integes and 0. It also
#outputs recurrences if their complexity is <=MaxC. Try:
#BookA2(2,10);
BookA2:=proc(M,MaxC) local NEG,POS,c1,i,gu1,gu2,gu11,gu21,t0,S:
t0:=time():
print(``):
print(`Enumerating Generalized Dyck paths, as well as the Sum of the Areas and Areas-Squared  with alphabets consisting of integers from`, -M, ` to `, M):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
c1:=0:
NEG:={seq(-i,i=1..M)}:
POS:={seq(i,i=1..M)}:
gu1:=powerset(NEG) minus {{}}:
gu2:=powerset(POS) minus {{}}:
print(``):
print(`-------------------------------------------------------------`):
print(``):

for gu11 in gu1 do
for gu21 in gu2 do
 S:=gu11 union gu21:

if igcd(op(S))=1 then

c1:=c1+1:

Paper2Ac(S,MaxC,c1):
print(``):
print(`-------------------------------------------------------------`):
print(``):
S:=S union {0}:
c1:=c1+1:
Paper2Ac(S,MaxC,c1):
fi:
od:
od:

print(``):
print(`--------------------------------`):
print(``):
print(`This ends this book that took`, time()-t0, `seconds to create `):
print(``):
end:


#ModeP(P,q): inputs a polynomial P in q and outputs the pair [degree,mode] where degree is its degree and mode is the power with the largest coefficient
#assuming that it is unimodal, it returns FAIL otherwise. Try:
#ModeP(1+3*q+2*q^2,q);
#ModeP(qSeqS({1,0,-1},q,10)[11],q);
ModeP:=proc(P,q)  local i,d,m:
d:=degree(P,q):
for i from 0 while coeff(P,q,i)<=coeff(P,q,i+1) do od:
m:=i:
for i from m to d-1 do
 if coeff(P,q,i)<coeff(P,q,i+1) then
  RETURN([d,FAIL]):
 fi:
od:
[d,m]:
end:


#end from qGDW.txt

#start Strict Walks

#



#EqGFtS(S,X,t): the algebraic equation satisfied by the generating function for the sequence enumerating STRICT generalized Dyck paths with set of steps S of length n not keeping track of the original steps. Try:
#EqGFtS({1,-1},X,t);
EqGFtS:=proc(S,X,t) local gu,x,f,g,eq,var,s1,eq1:
option remember:
gu:=MakeSys(f,g,x,S):
eq:=gu[1]:
var:=gu[2] minus {g[0,0]}:

eq:=subs({seq(x[s1]=t,s1 in S)},eq):


eq:={seq(op(1,eq1)-op(2,eq1),eq1 in eq)}:

eq:=subs(g[0,0]=X,Groebner[Basis](eq,plex(op(var),g[0,0]))[1]):
eq:=factor(eq):

eq:

end:


#qEqGFtS(S,X,t): Inputs a set of positive and negative integers S, outputs an algebraic equation for X=X(t) the
#generating function for the sun-of-the-areas under all generalized STRICT Dyck paths with steps {[1,s], s in S}. For the
#classical Dyck case do:
#qEqGFtS({1,-1},X,t);

qEqGFtS:=proc(S,X,t) local q,f,g,sys, eq, var,eq1,v,EQ0,EQ1,EQ2,EQ,VAR,i,Z,lu,aj,S1,s:
aj:=igcd(op(S)):

if aj>1 then
 S1:={seq(s/aj,s in S)}:
 RETURN(numer(subs(X=X/aj,qEqGFt(S1,X,t)))):
fi:
	sys:=qMakeSyst(f,g,t,q,S):
	eq:=sys[1]:
	var:=sys[2]: 

	eq1:=expand(subs( {seq(v = op(0, op(0, v))[op(1 .. 2, op(0, v)), 0](op(1, v)) + (q - 1)*op(0, op(0, v))[op(1 .. 2, op(0, v)), 1](op(1, v)), v in var)},eq) 
	): 
	var:=subs(q=1,var):  
	VAR:={seq(op(0,op(0,v))[op(1..2,op(0,v)),0](op(1,v)),v in var),seq(op(0,op(0,v))[op(1..2,op(0,v)),1](op(1,v)),v in var),seq(diff(op(0,op(0,v))[op(1..2,op(0,v)),0](op(1,v)),op(1,v)),v in var)}:  

	EQ0:=subs(q=1,eq1):
	EQ1:=subs(q=1,diff(eq1,q)):
	EQ2:=subs(q=1,diff(eq1,t)):
	EQ:=EQ0 union EQ1 union EQ2: 
	EQ:=subs({seq(D(op(0,op(0,v))[op(1..2,op(0,v)),0]) (op(1,v))=diff(op(0,op(0,v))[op(1..2,op(0,v)),0](op(1,v)),op(1,v)),v in var)},EQ):
	VAR:=[op(VAR minus {g[0,0,1](t)}),g[0,0,1](t)]:
	EQ:=subs({seq(VAR[i]=Z[i],i=1..nops(VAR))},EQ): 
       
	VAR:=subs({seq(VAR[i]=Z[i],i=1..nops(VAR))},VAR): 
	lu:=subs(Z[nops(VAR)] = X, Groebner[Basis](EQ, plex(op(VAR)))[1]):
	lu:=add(factor(coeff(lu,X,i))*X^i,i=0..degree(lu,X)):
     lu:



 RETURN(factor(lu)):
end: