######################################################################
## 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 Doron Zeilberger, Rutgers University ,                 #
##  zeilberg at math dot rutgers dot edu.                            # 
######################################################################

with(combinat):

print(`First Written: Oct. , 2022: 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(` Experimenting with the Area Statistics of  Generalized Dyck Paths`):

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



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(` AveSD, DecoW, DerEqP, DerEqPk, DerF, EqGFe, Findrec, GFpr, GFstx, IsStrict, MakeEqF,   MakeSys,  MakeSysE, NuWabn, PWabn, PWE, 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, EqP, SeqFromEq, SeqS, SeqSe `):

print(``):


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])=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])=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, followed by the factor that you have`):
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, outouts 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])=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])=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])=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])=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])=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])=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:

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

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

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:
#Paper({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, followed by the factor that you have
#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:
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: