#Anthony Zaleski 
#Recursively finding the generating function in Zeilberger's Feller paper:
#http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimPDF/feller.pdf



#read(`Moments3D.txt`):


#LookupProc(a): inputs procedure, outputs description
LookupProc:=proc(a)
if a=Fab then
	print(`Fab(a,b,t): inputs nonneg. integers a,b and variable`):
	print(`(or list of 4 variables) t. Uses recursion to find`):
	print(`sum_{walks from (0,0)->(a,b)} t[1]^a1 t[2]^a2 t[3]^a3 t[4]^a4,`);
	print(`where the ai's are the statistics from Feller.`):
elif a=FabLong then
	print(`FabLong(a,b,t): inputs nonneg. integers a,b and variable`):
	print(`(or list of 4 variables) t. Uses brute force to find`):
	print(`sum_{walks from (0,0)->(a,b)} t[1]^a1 t[2]^a2 t[3]^a3 t[4]^a4,`);
	print(`where the ai's are the statistics from Feller.`):
elif a=Help1 then
	print(`Help1(): the help function for secondary procedures.`):
elif a=Help then
	print(`Help(): the help function for primary procedures.`):
elif a=A then
	print(`A(L): inputs list of elements in {-1,1}, outputs the Feller statistics [a1,a2,a3,a4].`):
elif a=FabU then
	print(`FabU(a,b,t): same as Fab(a,b,t), but only counts walks ending with Up.`):
elif a=FabR then
	print(`FabR(a,b,t): same as Fab(a,b,t), but only counts walks ending with Right.`):
elif a=Fn then
	print(`Fn(n,t): sum_{walks of length n} t[1]^a1 t[2]^a2 t[3]^a3 t[4]^a4,`):
	print(`where the ai's are the statistics from Feller.`): 
elif a=ID then
	print(`ID(f,x,k): Given the generating function f`):
	print(`in a single variable x, outputs the list`):
	print(`[1st moment, variance, third moment about mean,...,kth moment about mean]`):
elif a=IDN then
	print(`IDN(f,x,k): Same as ID, but normalized.  E.g. try`):
	print(`evalf(IDN((1+x)^100,x,10))`):
elif a=CheckAnsatz then
	print(`CheckAnsatz(X,Y,f,var): inputs lists of n-tuples X and values Y, `):
	print(`and an expression f in the variables listed in var. Outputs true iff f interpolates X,Y`):
	print(`If n=1, it is OK to exclude around X-values and var:`):
	print(`E.g., CheckAnsatz([1,2,3],[1,4,9],x^2,x)=true.`):
elif a=GuessPoly1 then
	print(`GuessPoly1(X,Y,x,Deg): tries to find the UNIQUE degree Deg polynomial`):
	print(`in x that interpolates the lists the lists X,Y.  Returns FAIL if not possible.`):
elif a=GuessPoly then
	print(`GuessPoly(X,Y,x,MaxDeg:=10,MinDeg:=0): finds the polynomial in x with smallest degree between`):
	print(`MinDeg (0 by default) and MaxDeg (10 by default) which interpolates points given by`):
	print(`lists X and Y of x-values and y-values, resp.`):
	print(`If nops(X) is too small to make a reliable guess, returns FAIL.`):
elif a=IDs then
	print(`IDs(L,x,k,n,MaxDeg:=10,MinDeg:=0): given a list L=[f[1],...,f[N]] of generating functions`): 
	print(`in the variable x, outputs a list M=[M[1],...M[k]], where M[i] is an expression describing`):
	print(`the i^th moment of f[n] as a function of n.  (All moments except the first are about the mean.)`):
	print(`Each M[i] is a polynomial of degree between MinDeg (0 by default) and MaxDeg (10 by default).`):
	#print(`If it can't be done for all the moments up to the kth, goes up to the highest one possible.`):
elif a=IDsN then
	print(`IDsN(L,x,k,n,MaxDeg:=10,MinDeg:=0): normalized version of IDs.`):
elif a=IDsNL then
	print(`IDsNL(L,x,k,MaxDeg:=10,MinDeg:=0): IDs after normalizing and letting n->infinity.`):
elif a=CheckNormal then
	print(`CheckNormal(L): checks if the list L is of the form [0,1,0,3,0,15,...].`):
elif a=xD then
	print(`xD(f,x): the operator (f,x)->x*diff(f,x), used in computing moments.`):
elif a=ID2 then
	print(`ID2(f,x,y,k,stan:=true)): Inputs gen. fun. f in two variables x,y and positive integer k.`):
	print(`Outputs k+1 by k+1 list M s.t. M[i][j] is the i-1,j-1  standardized mixed moment of X,Y about their`):
	print(`means.  E.g., L[2][2] is the covariance.  set stan=false to not standardize.`):
	print(`E.g., try ID2(F3([10,10,10],[x,y,1$5]),x,y,1).`):
elif a=ID2s then
	print(`ID2s(L,x,y,k,n,MaxDeg:=10,MinDeg:=0): inputs list L=[f[1],..f[N]] of gen. fun's in x,y;`):
	print(`outputs k by k list M of expressions in n s.t. M[i][j] gives the i,j mixed moment`):
	print(`about the mean of f[n].`):
	print(`The expressions are polynomials with degree between MinDeg (0 by default) and MaxDeg (10 by default).`):
	#print(`If no such M exists, returns FAIL.`):
elif a=GuessPolyG then
	print(`GuessPolyG(X,Y,n,p,MaxDeg:=10,MinDeg:=0): Inputs lists X,Y of n-values and`):
	print(`y-values. Outputs a list P=[f_1,..,f_p] of polynomials in n  s.t. `):
	print(`Y(n)=f_i(n) if n=i mod p.`):
	print(`For example, GuessPolyG([1,2,3,4,5,6,7,8,9,10],[0,2,0,4,0,6,0,8,0,10],n,2)`):
	print(`returns [0,n].`):
elif a=IDsG then
	print(`IDsG(L,x,k,n,p,MaxDeg:=10,MinDeg:=0): given a list L=[f[1],...,f[N]] of generating functions `):
	print(`in the variable x, outputs a list M=[M[1],...M[k]], where M[i] is a list of p polynomials describing`):
	print(`the i^th moment of f[n] as a function of n.  (All moments except the first are about the mean.)`):
	print(`Each M[i] is a polynomial of degree between MinDeg (0 by default) and MaxDeg (10 by default).`):  
	print(`For example, try`):
	print(`IDsG([seq(Fn(n,[t,1,1,1]),n=1..20)],t,5,n,2);`):
elif a=IDsGN then
	print(`IDsGN(L,x,k,n,p,MaxDeg:=10,MinDeg:=0): normalized version of IDsG.`):
	print(`For example, try IDsGN([seq(Fn(n,[t,1,1,1]),n=1..20)],t,5,n,2);`):
elif a=IDsGNL then
	print(`IDsGNL(L,x,k,p,MaxDeg:=10,MinDeg:=0): normalized version of IDsG,`):
	print(`after letting n->infinity.`):
	print(`For example, try IDsGNL([seq(Fn(n,[t,1,1,1]),n=1..20)],t,5,2);`):
elif a=IsDecreasing then
	print(`IsDecreasing(L): checks if the list L decreases strictly.`):
	print(`For example, IsDecreasing([3,2,1]) returns true, but`):
	print(`IsDecreasing([3,2,2]) returns false.`):
elif a=F1 then
	print(`F1(L,t): the generating function for the d-dimensional generalization`):
	print(`of a_1, the number of losing times.  For example,`):
	print(`F1([n,n],0)`):
	print(`returns the nth Catalan number.`):
elif a=NballsDboxes then
	print(`NballsDboxes(n,d): set of all lists of d nonneg. integers`):
	print(`summing to n.  For example,`):
	print(`NballsDboxes(3,2) returns {[0,3], [1,2], [2,1], [3,0]}.`):
elif a=F1n then
	print(`F1n(n,t,d): the sum over d-dimensional walks of length n`):
	print(`of t^a_1, where a_1 is the generalized number of losing times.`):
	print(`For example,`):
	print(`F1n(4,t,2) = Fn(4,[t,1,1,1]).`):
elif a=dTerms then
	print(`dTerms(d,var): outputs set of products of powers of variables`):
	print(`in the set var having total degree exactly d.  E.g., `):
	print(`dTerms(2,{x,y})={x^2,y^2,xy}.`):
elif a=GuessRat1 then
	print(`GuessRat1(X,Y,d,var)inputs a list X of n-tuples and`):
	print(`a list Y of values at the n-tuples.`):
	print(`Outputs a guess for the rational function P/Q of`): 
	print(`degree d relating X and Y, as an expression in terms `):
	print(`of the elements of list var, where nops(var)=n.`):
	print(`E.g.,`):
	print(`GuessRat1([seq(seq([i,j],j=1..5),i=1..5)],[seq(seq((i+j)/i,j=1..5),i=1..5)],2,[x,y]);`):
	print(`returns (x+y)/x.`):
elif a=GuessRat then
	print(`GuessRat(X,Y,dmax,var): inputs a list X of n-tuples and`):
	print(`a list Y of values at the n-tuples.`):
	print(`Outputs a guess for the rational function P/Q of`): 
	print(`degree <=dmax relating X and Y, as an expression in terms `):
	print(`of the elements of list var, where nops(var)=n.`):
	print(`E.g.,`):
	print(`GuessRat([seq(seq([i,j],j=1..5),i=1..5)],[seq(seq((i+j)/i,j=1..5),i=1..5)],2,[x,y]);`):
	print(`returns (x+y)/x.`):
	print(`If nops(var)=1, OK to let X be a list of NUMBERS.`):
elif a=WinLose then
	print(`WinLose(a,b,x,y): the sum over walks to (a,b) of x^a*y^b, where`);
	print(`a is the number of winning times and b is the number of losing times.`):
	print(`E.g., try WinLose(2,3,x,y)`):
elif a=Kd then
	print(`Kd(d,a,b): the degree d coefficient in the generating function for the`):
	print(`number of losing times of a walk to (a,b).  E.g., Kd(0,i,i)`):
	print(`returns the ith Catalan number.`):
elif a=Cd then 
	print(`Cd(d,n,i): Kd(d,n,n+i)/binomial(2*n+i,n). E.g., Cd(0,i,0) returns 1/(i+1).`):
elif a=GuessCd then
	print(`GuessCd(d,n0,nInc,nSteps,i0,iInc,iSteps,n,i,dmax): for fixed d, guesses Cd(d,n,i) `):
	print(`as a rational function in terms of n and i. For input data, uses a grid in the `):
	print(`n-i plane specified by n0 (initial n), nInc (step size in n direction),`): 
	print(`nSteps (number of steps in n direction), and similarly for i0, iInc, iSteps.`):
	print(`Returns a rational function of degree <=dmax in variables n and i.`):
	print(`E.g., try GuessCd(0,0,1,10,0,1,10,n,i,10);`):
elif a=FactorCoef then
	print(`FactorCoef(p,x): factors the coefficients of polynomial p(x).`):
	print(`E.g., try FactorCoef(4*x^2+6,x)`):
elif a=F2 then
	print(`F2(L,t): sum over walks to L of t^(number of visits to x1>=x2>=...)`):
	print(`E.g., try F1([3,3,3],t).`):
elif a=Zinn then
	print(`Zinn(resh): Zinn-Justin's method to estimate`):
	print(`the [theta,mu] such that`):
	print(`resh[i] is appx. Const*mu^i*i^theta`):
	print(`For example, try:`):
	print(`Zinn([seq(5*i*2^i,i=1..30)]);`):
elif a=Moments3D then
	print(`Moments3D.txt: contains lists of the 1st through 5th moments of`):
	print(`F2([n,n,n],t),n=1..100.  E.g., try Mom1;`):
elif a=FitAnsatz then
	print(`FitAnsatz(X,Y,x,ansatz,var): inputs lists X and Y of length`):
	print(`at least nops(var), symbol x, ansatz in terms of x and the variables`):
	print(`in the set var, and outputs the ansatz which EXACTLY interpolates`):
	print(`the first nops(var) data points.`):
	print(`E.g., try FitAnsatz([1,2,3],[2,3,4],x,a*x+b,{a,b}).`):
elif a=LogLogLinear then
	print(`LogLogLinear(X,Y,x,XY:=0): finds least squares linear fit to log(X), log(Y).`):
	print(`If XY=0, outputs the fit as an expression ax+b.`):
	print(`If XY=1, outputs [fit,log(X),log(Y)].`):
	print(`E.g., try LogLogLinear([1,2,3],[1,4,9],x,1);`):
elif a=F3 then
	print(`F3(L,T): Inputs lists L=[x,y,z] and T=[t123,t132,t213,t231,t312,t321,t0]`):
	print(`Outputs generating function for walks to L where t123 keeps track of`):
	print(`visits to x>y>z, etc., and t0 keeps track of other visits.`):
	print(`E.g., try F3([10,10,10],[x,y,1$5]).`):
elif a=F2G then
	print(`F2G(L,t,S): inputs list L=[a,b], variable t, and set of steps S.`):
	print(`Outputs generalized Feller polynomial.  E.g., try`):
	print(`F2G([3,3],t,{[1,0],[0,1]})`):
elif a=F3G then
	print(`F3G(L,t,S): inputs list L=[a,b,c], variable (list) t, and set of steps S.`):
	print(`Outputs generalized Feller polynomial in t[1]..t[7].  E.g., try`):
	print(`F3G([3,3,3],t,{[1,0,0],[0,1,0],[0,0,1]})`):
elif a=GuessPower then
	print(`GuessPower(X,Y,Round:=true): inputs lists X,Y, outputs guess p s.t.`):
	print(`Y[i] is appx. C*X[i]^p.  If Round=true, rounds p.`):
	print(`E.g., try GuessPower([1,2,3],[1,4,10]);`): 
elif a=Aitken1 then
	print(`Aitken1(a1,a2,a3): Returns one number, the Aitken transformation.`):
elif a=FellerBook2D then
	print(`FellerBook2D(S,M,K1,K2): inputs a set of steps S, natural numbers M,K1<K2.`):  
	print(`Prints a book about moments 1-M of 2D walks from (0,0)->(n,n) with steps from S.`):
	print(`Uses n=K1..K2 to guess moment behavior. `):
	print(`If verbose is 1 or true, outputs book; otherwise, outputs list with`):
	print(`entries of form [set, c_1, c_2,...,c_M].`):
	print(`E.g., try FellerBook2D({[0,1],[1,0],[1,1]},3,10,20);`);
elif a=FellerBook2DSafe then 
	print(`FellerBook2DSafe(S,M,K1,K2,K3,K4,verbose:=1): same as FellerBook2D,`):
	print(`but guesses the constants using N=K1..K2 and N=K3..K4 and only keeps`):
	print(`the digits that agree.`):
elif a=ChungFellerGF then
	print(`ChungFellerGF(t,k,pgf:=0,center:=0): The weight enumerator of`):
	print(`t^(#losing times) for walks from (0,0)->(k,k).`):
	print(`Both t and k may be numbers or symbols.`):
	print(`pgf=1 returns the PGF; center=0 divides by t^k.`):
	print(`E.g., try ChungFellerGF(t,3);`):
elif a=ChungFellerMoment then
	print(`ChungFellerMoment(m,k,standard:=1,centered:=1): the mth moment of the`):
	print(`C-F generating function for walks to (k,k) (k may be a symbol).`):
	print(`Options standard and centered specify if the moment is standardized/centered.`):
	print(`E.g., ChungFellerMoment(2,k,0,1); computes variance.`):
elif a=ForwardKingGF then
	print(`ForwardKingGF(z,t): the sum over walks from the origin to y=x`):
	print(`with step set {[1,0],[0,1],[1,1]} of z^k t^a, where the walk`):
	print(`ends at (k,k) and is losing a times.`):
	print(`E.g., try ForwardKingGF(z,t);`):
else
	print(`No such procedure exists.`):
fi:
end:




Help:=proc():
if nargs=0 then
	print(`Feller.txt contains the following primary procedures:`):
	print(`EXPERIMENTAL:`):
	print(`Fab,Fn,F1,F1n,F2,F3,F2G,F3G,FellerBook2D`):
	print(`ANALYTIC:`):
	print(`ChungFellerGF,ChungFellerMoment,ForwardKingGF`):
	print(`MOMENTS:`):
	print(`ID,IDN,IDs,IDsN,IDsNL,ID2,ID2s,Moments3D,IDsG,IDsGN,IDsGNL.`):	
	print(`---`):
	print(`Help1 lists secondary procedures.`):
else
	LookupProc(args[1]):
	
fi:
end:


Help1:=proc():
if nargs=0 then
	print(`Feller.txt contains the following secondary procedures (Help() lists the primary procedures):`):
	print(`A,FabLong,FabR,FabU`):
	print(`CheckAnsatz,GuessPoly1,GuessPoly,GuessPolyG,xD`):
	print(`NballsDboxes `):
	print(`IsDecreasing `):
	print(`dTerms,GuessRat1, GuessRat`):
	print(`Zinn,FitAnsatz,LogLogLinear,GuessPower,Aitken1`):
	print(`Kd, Cd`):
	print(`FactorCoef `):
else
	LookupProc(args[1]):
	
fi:
end:




with(combinat):


######################################
######################################



A:=proc(L) local i,S,a1,a2,a3,a4:

S:=0: a1:=0: a2:=0: a3:=0: a4:=0:

for i from 1 to nops(L) do
	S:=S+L[i]:
	if S<0 or (S=0 and L[i]=1) then
		a1:=a1+1:
	fi:

	if S=0 then
		a2:=a2+1:
		a3:=i:
		if i<nops(L) then
			if L[i]=L[i+1] then
				a4:=a4+1:
			fi:
		fi:
	fi:

	if S<0 then
		a3:=0:
	fi:

od:
[a1,a2,a3,a4]:
end:



Fab4Long:=proc(a,b,t) local W,i,f,L:
W:=permute([-1$a,1$b]):
f:=0:
for i from 1 to nops(W) do
	L:=A4(W[i]):
	f:=f+t[1]^L[1]*t[2]^L[2]*t[3]^L[3]*t[4]^L[4]:
od:
f:
end:


FabR:=proc(a,b,t)
local T:
option remember:

T:=[t[1],t[2],1,t[4]]:

if a=0 then
	return 0:
fi:

if b=0 then
	return t[1]^a:
fi:

if a<0 or b<0 then
	return 0:
fi:

if b=a then
	return expand( t[2]*t[3]^(a+b)*(Fab(a-1,b,T)) ):
fi:

if b>a then 
	return Fab(a-1,b,t):
fi:

if b=a-1 then
	return expand( t[1]*t[4]*FabR(a-1,b,T)+t[1]*FabU(a-1,b,T) ):
fi:

expand( t[1]*Fab(a-1,b,T) ):
end:



FabU:=proc(a,b,t)
local T:
option remember:

T:=[t[1],t[2],1,t[4]]:

if b=0 then
	return 0:
fi:

if a=0 then
	return 1:
fi:

if a<0 or b<0 then
	return 0:
fi:

if b=a then
	return expand( t[1]*t[2]*t[3]^(a+b)*(Fab(a,b-1,T)) ):
fi:

if b=a+1 then 
	return expand( FabR(a,b-1,t) + t[4]*FabU(a,b-1,t) ):
fi:

if b>a+1 then
	return Fab(a,b-1,t):
fi:

expand( t[1]*Fab(a,b-1,T) ):

end:



Fab:=proc(a,b,t)
option remember:
if a=0 and b=0 then
	return 1:
fi:
FabR(a,b,t)+FabU(a,b,t):
end:


Fn:=proc(n,t) local i:
add(Fab(i,n-i,t),i=0..n):
end:







#Finding moments symbolically 
######################################
######################################


CheckAnsatz:=proc(X,Y,f,var) local neq,j,nvar,i:

neq:=nops(X):
nvar:=nops(var):

if not type(X[1],`list`) then
	return CheckAnsatz([seq([X[i]],i=1..neq)],Y,f,[var]):
fi:

for i from 1 to neq do
	if Y[i]<>subs({seq(var[j]=X[i][j],j=1..nvar)},f) then
		return false:
	fi:
od:
true:
end:


GuessPoly1:=proc(X,Y,x,Deg) local eq,var,a,P,i:
P:=add(x^i*a[i],i=0..Deg):
var:={seq(a[i],i=0..Deg)}:
eq:={seq(subs(x=X[i],P)=Y[i], i=1..Deg+1)}:

var:=solve(eq,var):
P:=subs(var,P):

if var=NULL or CheckAnsatz(X,Y,P,x)=false then 
	return FAIL: 
fi:
P:
end:


GuessPoly:=proc(X,Y,x,MaxDeg:=10,MinDeg:=0) local d,P:
if nops(X)<>nops(Y) then
	print(`ERROR in GuessPoly: X,Y must be lists of the same length.`):
	return(FAIL):
fi:

for d from MinDeg to MaxDeg do
	if nops(X)<d+4 then
		#print(`ERROR in GuessPoly: the lists X and Y are too short`):
		#print(`to make a reliable guess.  Try more data points.`):
		return FAIL:
	fi:

	P:=GuessPoly1(X,Y,x,d):
	if P<>FAIL then
		return P:
	fi:
od:
FAIL:
end:



GuessPolyG:=proc(X,Y,n,p,MaxDeg:=10,MinDeg:=0) local X1,Y1,i,j,P:
if nops(X)<>nops(Y) or nops(X)<p then
	print(`ERROR in GuessPolyG: need X,Y to specify the`):
	print(`x and y values of at least p data points.`):
	return FAIL:
fi:

P:=[]:

for i from 1 to p do
	X1:=[]: Y1:=[]:
	for j from i by p to nops(X) do
		X1:=[op(X1),X[j]]:
		Y1:=[op(Y1),Y[j]]:
	od:
	P:=[op(P), GuessPoly(X1,Y1,n,MaxDeg,MinDeg)]:
od:
P:
end:



#Credit: http://www.math.rutgers.edu/~zeilberg/EM15/C6.txt

ID:=proc(e,x,k) local i,av,M, f:

f:=e/subs(x=1,e):

av:=subs(x=1,diff(f,x)):

M:=[av]:


f:=f/x^av:


f:=x*diff(f,x):
for i from 2 to k do
f:=x*diff(f,x):
M:=[op(M),subs(x=1,f)]:
od:
 
M:
end:


IDN:=proc(e,x,k) local i,av,M, f:

if e=0 then return([0$k]): fi:

f:=e/subs(x=1,e):


av:=subs(x=1,diff(f,x)):

M:=[av]:


f:=f/x^av:


f:=x*diff(f,x):
for i from 2 to k do
f:=x*diff(f,x):
M:=[op(M),subs(x=1,f)]:
od:

if M[2]<>0 then 
 return [M[1],M[2],seq(M[i]/M[2]^(i/2),i=3..k)]:
fi:

[seq(M[i],i=1..k)]:
end:




IDs:=proc(L,x,k,n,MaxDeg:=10,MinDeg:=0) local Mn,N,M,i,j,X,Y,P:
N:=nops(L):
Mn:=[seq(ID(L[j],x,k),j=1..N)]:
M:=[]:
X:=[seq(j,j=1..N)]:

for i from 1 to k do
	Y:=[seq(Mn[j][i], j=1..N)]:
	P:=GuessPoly(X,Y,n,MaxDeg,MinDeg):

	#if P=FAIL then return M: fi:

	M:=[op(M),P]:
od:
if FAIL in M then
	print(`ERROR in IDs: the list L was not long enough`):
	print(`to reliably guess all the moments.`):
fi:
M:
end:



IDsG:=proc(L,x,k,n,p,MaxDeg:=10,MinDeg:=0) local Mn,N,M,i,j,X,Y,P:
N:=nops(L):
Mn:=[seq(ID(L[j],x,k),j=1..N)]:
M:=[]:
X:=[seq(j,j=1..N)]:

for i from 1 to k do
	Y:=[seq(Mn[j][i], j=1..N)]:
	P:=GuessPolyG(X,Y,n,p,MaxDeg,MinDeg):

	#if P=FAIL then return M: fi:

	M:=[op(M),P]:
od:
M:
end:



IDsGN:=proc(L,x,k,n,p,MaxDeg:=10,MinDeg:=0) local M,i,sig:
M:=IDsG(L,x,k,n,p,MaxDeg,MinDeg):
for i from 1 to p do
	sig[i]:=sqrt(M[2][i]):
od:
[0,1,seq([seq(normal(M[j][i]/sig[i]^j),i=1..p)],j=3..k)]:
end:



IDsGNL:=proc(L,x,k,p,MaxDeg:=10,MinDeg:=0) local M,i,sig,n:
M:=IDsG(L,x,k,n,p,MaxDeg,MinDeg):
for i from 1 to p do
	sig[i]:=sqrt(M[2][i]):
od:
[0,1,seq([seq(limit(normal(M[j][i]/sig[i]^j),n=infinity),i=1..p)],j=3..k)]:
end:



IDsN:=proc(L,x,k,n,MaxDeg:=10,MinDeg:=0) local M,sig:
M:=IDs(L,x,k,n,MaxDeg,MinDeg):

if M=FAIL then return FAIL: fi:

sig:=sqrt(M[2]):
[0,seq(normal( M[i]/sig^i ),i=2..nops(M) )]:
end:


IDsNL:=proc(L,x,k,MaxDeg:=10,MinDeg:=0) local M,sig,n:
M:=IDs(L,x,k,n,MaxDeg,MinDeg):

if M=FAIL then return FAIL fi:

sig:=sqrt(M[2]):
[0,seq(limit( normal(M[i]/sig^i) ,n=infinity ),i=2..nops(M) )]:
end:



CheckNormal:=proc(L) local i,j:
for i from 1 by 2 to nops(L) do
	if L[i]<>0 then
		return false:
	fi:
od:

for i from 2 by 2 to nops(L) do
	if L[i]<>mul(2*j+1,j=1..i/2-1) then
		return false:
	fi:
od:
true:
end:


xD:=(f,x)->x*diff(f,x):

ID2:=proc(f,x,y,k,stan:=true) local s,g,L,i,j,Ex,Ey,gx,gxy,MomX,MomY,vX,vY:
s:={x=1,y=1}:
g:=f/subs(s,f):

MomX:=ID(subs(y=1,g),x,k):
MomY:=ID(subs(x=1,g),y,k):

Ex:=MomX[1]:
Ey:=MomY[1]:
if k>1 then
	vX:=MomX[2]:
	vY:=MomY[2]:
fi:

L[0][0]:=1:
for i from 1 to k do
	if (stan=true or stan=1) and i>=3 then
		L[i][0]:=MomX[i]/vX^(i/2):
		L[0][i]:=MomY[i]/vY^(i/2):
	else
		L[i][0]:=MomX[i]:
		L[0][i]:=MomY[i]:
	fi:
od:


#Ex:=subs(s,xD(g,x)):
#Ey:=subs(s,xD(g,y)):


g:=g/x^Ex/y^Ey:
gx:=g:

for i from 1 to k do
	gx:=xD(gx,x):
	gxy:=gx:
	for j from 1 to k do
		gxy:=xD(gxy,y):
		if (stan=true or stan=1)  then
			L[i][j]:=subs(s,gxy)/vX^(i/2)/vY^(j/2):
		else
			L[i][j]:=subs(s,gxy):
		fi:
		
	od:
od:

evalf( [seq([seq(L[i][j],j=0..k)], i=0..k)] ):
end:



ID2s:=proc(L,x,y,k,n,MaxDeg:=10,MinDeg:=0) local i,j,m,Mn,N,M,X,Y,flag:
N:=nops(L):
Mn:=[seq(ID2(L[m],x,y,k) , m=1..N)]:

X:=[seq(m, m=1..N)]:
flag:=false:
for i from 1 to k do
	for j from 1 to k do
		Y:=[seq(Mn[m][i][j], m=1..N)]:
		M[i][j]:=GuessPoly(X,Y,n,MaxDeg,MinDeg):

		if M[i][j]=FAIL then flag:=true: fi:
	od:
od:
if flag then 
	print(`ERROR in ID2s: the list L was not long enough`):
	print(`to reliably guess all the moments.`):
fi:
[seq([seq(M[i][j],j=1..k)],i=1..k)]:
end:



######################################
#Higher-dimensional walks
######################################




IsDecreasing:=proc(L) local i:
for i from 1 to nops(L)-1 do
	if L[i+1]>=L[i] then
		return false:
	fi:
od:
true:
end:




F1:=proc(L,t) local d,f,L1,i,j:
option remember:

d:=nops(L):

if `or`(seq(L[i]<0, i=1..d)) then return 0: fi:

if L=[0$d] then return 1: fi:

f:=0:

for i from 1 to d do 
	L1:=[seq(L[j],j=1..i-1), L[i]-1, seq(L[j],j=i+1..d)]:
	if  IsDecreasing(L1) or IsDecreasing(L) then 
	    #IsDecreasing(L) or ( nops({op(L)})=1 and IsDecreasing(L1)) then 
	    #IsDecreasing(L) then	    
		f:=f+F1(L1,t)*t:
	else
		f:=f+F1(L1,t):
	fi:
od:
expand(f):
end:


NballsDboxes:=proc(n,d) local S,L,P,i,j:
option remember:
if n=0 then return { [0$d] }: fi:

S:=NballsDboxes(n-1,d):
P:={}:

for L in S do
	for i from 1 to d do
		P:=P union { [seq(L[j],j=1..i-1), L[i]+1, seq(L[j],j=i+1..d)] }
	od:
od:
P:
end:


F1n:=proc(n,t,d) local S,L:
option remember:
S:=NballsDboxes(n,d):
add(F1(L,t), L in S):
end:



######## Using different conventions

IsNoninc:=proc(L) local i:
for i from 1 to nops(L)-1 do
	if L[i+1]>L[i] then
		return false:
	fi:
od:
true:
end:

F2:=proc(L,t) local d,f,L1,i,j:
option remember:

d:=nops(L):

if `or`(seq(L[i]<0, i=1..d)) then return 0: fi:

if L=[0$d] then return 1: fi:

f:=0:


for i from 1 to d do 
	L1:=[seq(L[j],j=1..i-1), L[i]-1, seq(L[j],j=i+1..d)]:
	f:=f+F2(L1,t):
od:

if IsNoninc(L) then
	f:=t*f:
fi:
expand(f):
end:



#######################################
#Guessing multivariable rational functions!
########################################

dTerms:=proc(d,var,degBounds:=0) local a,i,d1,var1,x,S,s,S1: option remember:
if d=0 then
	return {1}:
fi:

S:=`union`(  seq({op(expand(a*[op(dTerms(d-1,var))]))}, a in var)  ):

if degBounds<>0 then
	if (not type(degBounds,`list`)) or 
	    nops(degBounds)<>nops(var) then
		print(`dTerms ERROR: degBounds must be`):
		print(`0 or list of length nops(var).`):
		return FAIL:
	fi:
	S1:={}:
	for s in S do
		if `and`( seq(degree(s,var[i])<=degBounds[i],i=1..nops(var)) ) then
			S1:=S1 union {s}:
		fi:
	od:

	return S1:
fi:

S:
end:



#Note: the required number of equations
#increases in both d and the numbers of vars
#for |var|=1 we seem to need 2*d+6 equations.

GuessRat1:=proc(X,Y,d,var,degBounds:=0,restrict:=false) local S,i,T,a,b,P,Q,t,eq,coeff,j,n,f,neq,P0,Q0, S0:
T:=`union`(seq(dTerms(i,[op(var)],degBounds),i=0..d)): 
n:=nops(T):
P:=add(a[i]*T[i], i=1..n):
Q:=add(b[i]*T[i], i=1..n):
coeff:={ seq(a[i], i=1..n), seq(b[i], i=1..n) }:

if restrict then
	for i from 1 to n do
		if subs(var[2]=0,T[i])<>0 then
			Q:=subs(b[i]=a[i],Q):
			coeff:=coeff minus {b[i]}:
		fi:
	od:
fi:

if nops(var)=1 then
	neq:=2*d+6:
else
	neq:=nops(Y): #nops(coeff)+6
fi:

if nops(Y) < neq then
	print(`GuessRat1 error: not enough data points.`):
	return FAIL,neq:
fi:

eq:={add(b[i],i=1..n)=1}: # to avoid multiple solutions


for i from 1 to neq do
	S:={seq(var[j]=X[i][j],j=1..nops(X[1]))}:
	eq:=eq union {subs(S,P)-Y[i]*subs(S,Q)}:
od:




S:=solve(eq,coeff):
if S=NULL then
	return FAIL:
fi:

f:=normal(subs(S,P/Q)):

if neq<nops(Y) and CheckAnsatz(X,Y,f,var)=false then
	return FAIL:
fi:

f:
end:



GuessRat:=proc(X,Y,dmax,var,degBounds:=0,restrict:=false) local d,cur,neq,nvar:

neq:=nops(Y):
nvar:=nops(var):

if not type(X[1],`list`) then
	return GuessRat([seq([X[i]],i=1..neq)],Y,dmax,[op(var)],degBounds,restrict):
fi:

cur:=GuessRat1(X,Y,0,var,degBounds,restrict):

for d from 1 to dmax while nops([cur])=1 and cur=FAIL do
	cur:=GuessRat1(X,Y,d,var,degBounds,restrict):
od:
cur:
end:


C:=proc(n) option remember:
(2*n)!/(n+1)!/n!:
end:




#####################
#Symmetric GF for losing and winning times

WinLose:=proc(a,b,x,y) 
option remember:

if a<0 or b<0 then return 0: fi:

if a=0 and b=0 then return 1: fi:

if b=a then return expand(x*WinLose(a-1,b,x,y)+y*WinLose(a,b-1,x,y)): fi:

if b<a then return expand(y*(WinLose(a-1,b,x,y)+WinLose(a,b-1,x,y))): fi:

expand(x*(WinLose(a-1,b,x,y)+WinLose(a,b-1,x,y))):

end:



############################
#Guessing coefficients of the GF of the number of losing times
############################

Kd:=proc(d,a,b) local t: option remember:
coeff(Fab(a,b,[t,1,1,1]),t,d):
end:

Cd:=proc(d,n,i) option remember: local t:
coeff(Fab(n,n+i,[t,1,1,1])/Fab(n,n+i,[1,1,1,1]),t,d):
end:


GuessCd:=proc(d,n0,nInc,nSteps,i0,iInc,iSteps,n,i,dmax,degBounds:=0,restrict:=false) 
local X,Y,nn,ii:
X:=[seq(seq([n0+nn*nInc, i0+ii*iInc],nn=0..nSteps),ii=0..iSteps)]:
Y:=[seq(seq(Cd(d,n0+nn*nInc, i0+ii*iInc),nn=0..nSteps),ii=0..iSteps)]:
GuessRat(X,Y,dmax,[n,i]):
end:



FactorCoef:=proc(p,x) local i:
add(x^i*ifactor(coeff(p,x,i)),i=0..degree(p,x)):
end:



##########################
#Zinn from http://www.math.rutgers.edu/~zeilberg/tokhniot/Zinn


sn:=proc(resh,n1):
-1/log(op(n1+1,resh)*op(n1-1,resh)/op(n1,resh)^2):
end:
 
#Zinn(resh): Zinn-Justin's method to estimate
#the [theta,mu] such that
#resh[i] is appx. Const*mu^i*i^theta
#For example, try:
#Zinn([seq(5*i*2^i,i=1..30)]);
Zinn:=proc(resh)
local s1,s2,theta,mu,n1,i:
if nops({seq(sign(resh[i]),i=1..nops(resh))})<>1 then
 RETURN(FAIL):
fi:

n1:=nops(resh)-1:
s1:=sn(resh,n1):
s2:=sn(resh,n1-1):
theta:=evalf(2*(s1+s2)/(s1-s2)^2):
mu:=evalf(sqrt(op(n1+1,resh)/op(n1-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
[theta,mu]:
end:

##########################
#Asymptotics
##########################

with(CurveFitting):

FitAnsatz:=proc(X,Y,x,ansatz,var) local eq:
if nops(X)<>nops(Y) or nops(X)<nops(var) then
	print(`ERROR in FitAnsatz: need nops(X)=nops(Y)>=nops(var).`):
fi:
eq:={seq(subs(x=X[i],ansatz)=Y[i],i=1..nops(var))}:
subs(solve(eq,var),ansatz):
end:

LogLogLinear:=proc(X,Y,x,XY:=0) local XX,YY,i,a,b:
if nops(X)<>nops(Y) then
	print(`ERROR in LogLogLinear: X,Y must be lists of same length.`):
fi:
XX:=[seq(evalf(log(X[i])),i=1..nops(X))]:
YY:=[seq(evalf(log(Y[i])),i=1..nops(Y))]:
if XY=1 then
	return [LeastSquares(XX,YY,x,`curve`=a*x+b), XX, YY]:
fi:
LeastSquares(XX,YY,x,`curve`=a*x+b):
end:

GuessPower:=proc(X,Y,Round:=true) local XX,YY,i,x,p:
XX:=[]: YY:=[]:
for i from 1 to nops(X) do
	if X[i]>0 and Y[i]>0 then
		XX:=[op(XX),X[i]]:
		YY:=[op(YY),Y[i]]:
	fi:
od:
if nops(XX)=0 then return 0 fi:

p:=coeff(LogLogLinear(XX,YY,x),x,1):

if Round=true or Round=1 then
	return round(p):
fi:
p:
end:

Aitken1:=proc(a1,a2,a3) local Denom:
Denom:=(a3-2*a2+a1):
if Denom<>0 then
	return(a3-(a3-a2)^2/Denom):
fi:
a3:
end:

#######################

#T=[t123,t132,t213,t231,t312,t321,t0]

F3:=proc(L,T) local f,i,j,L1:
option remember:
if L[1]<0 or L[2]<0 or L[3]<0 then return 0: fi:

if L=[0,0,0] then return 1: fi:

f:=0:

for i from 1 to 3 do 
	L1:=[seq(L[j],j=1..i-1), L[i]-1, seq(L[j],j=i+1..3)]:
	if IsDecreasing([L1[1],L1[2],L1[3]]) then
		f:=f+F3(L1,T)*T[1]:
	elif IsDecreasing([L1[1],L1[3],L1[2]]) then
		f:=f+F3(L1,T)*T[2]:
	elif IsDecreasing([L1[2],L1[1],L1[3]]) then
		f:=f+F3(L1,T)*T[3]:
	elif IsDecreasing([L1[2],L1[3],L1[1]]) then
		f:=f+F3(L1,T)*T[4]:
	elif IsDecreasing([L1[3],L1[1],L1[2]]) then
		f:=f+F3(L1,T)*T[5]:
	elif IsDecreasing([L1[3],L1[2],L1[1]]) then
		f:=f+F3(L1,T)*T[6]:
	else
		f:=f+F3(L1,T)*T[7]:
	fi:
od:
expand(f):
end:

##################################
#General step allowances
################################## 


IDRegion:=proc(L)
if IsDecreasing([L[1],L[2],L[3]]) then
	return 1:
elif IsDecreasing([L[1],L[3],L[2]]) then
	return 2:
elif IsDecreasing([L[2],L[1],L[3]]) then
	return 3:
elif IsDecreasing([L[2],L[3],L[1]]) then
	return 4:
elif IsDecreasing([L[3],L[1],L[2]]) then
	return 5:
elif IsDecreasing([L[3],L[2],L[1]]) then
	return 6:
else
	return 7:
fi:
end:


F2G:=proc(L,t,S) local IsUnder, s, L1,f,a:
option remember:

if L=[0,0] then return 1: fi:
if L[1]<0 or L[2]<0 then return 0 fi:

IsUnder:=evalb(L[1]>L[2]):
f:=0:
for s in S do
	L1:=L-s:
	if IsUnder or L1[1]>L1[2] then
		a:=t:
	else
		a:=1:
	fi:
	f:=f+a*F2G(L1,t,S):
od:
expand(f):
end:



F3G:=proc(L,t,S) local s,f,L1:
option remember:
if L=[0$3] then return 1: fi:
if L[1]<0 or L[2]<0 or L[3]<0 then return 0: fi:

f:=0:

for s in S do
	L1:=L-s:
	f:=f+t[IDRegion(L1)]*F3G(L1,t,S):
od:
expand(f):
end:

#############################
#Feller Books
############################

PowerAnsatz:=proc(m) 
if m=1 then return 1:
elif m=2 then return 2:
else return 0: fi:
end:




FellerBook2D:=proc(S,M,K1,K2,verbose:=1) local P,p,i,j,Div,Moms,t,k,m,X,Y,N,StartTime,Pow,C,
	DataList, Cur:

StartTime:=time():

if verbose=true or verbose=1 then
	Div:=`------------------------------------`:

	print(`WALKS IN THE PLANE FROM (0,0) TO (N,N)`): 
	print(``):
	print(Div):
	print(`This paper is concerned with walks from (0,0) to (N,N) where`):
	print(`steps from some subset of`,S,`are allowed.`):
	print(``):
	print(`Given a nonempty subset E of`,S,`define the generating function`):
	print(`f(N,E;t):=sum t^a(w), where w ranges over walks from (0,0)->(N,N)`):
	print(`with steps taken from E, and a is the statistic which counts`):
	print(`the number of times a walk visits {y<x} or visits {y=x} coming from {y<x}.`):
	print(``):
	print(`For each E, we generate f(N,E;t) for N=`,K1,K1+1,`...`,K2):
	print(`and use this data to conjecture how moments 1 through`,M):
	print(`behave as functions of N.`):
	print(Div):
	print(``):
else
	DataList:=[]:
fi:

P:=choose(S) minus {{}}:
i:=1:
X:=[seq(j,j=K1..K2)]:
N:=nops(X):

for p in P do
	Moms:=[seq(evalf(IDN(evalf( F2G([k,k],t,p) ),t,M)),k=K1..K2)]:
	for m from 1 to M do
		Y:=[seq(Moms[j][m],j=1..N)]:
		Pow:=PowerAnsatz(m):#GuessPower(X,Y):
		C[m]:=Aitken1(seq(Y[j]/X[j]^Pow,j=N-2..N)):
	od:
	if {seq(C[j],j=1..M)}<>{0.} then
		if verbose=1 or verbose=true then
			print(`Theorem`,i,`: `):
			print(`When steps from`,p,`are allowed:`):
		else
			Cur:=[p]:
		fi:
			for m from 1 to M do
				Pow:=PowerAnsatz(m):
				if verbose=1 or verbose=true then
					if Pow=0 then
						print(`Moment`,m,`behaves like`,C[m]):
					elif Pow=1 then
						print(`Moment`,m,`behaves like`,C[m],`n`):
					else
						print(`Moment`,m,`behaves like`,C[m],`n^`,Pow):
					fi:
				else
					Cur:=[op(Cur),C[m]]:
				fi:
		od:
		if verbose=1 or verbose=true then
			print(Div):
		else
			DataList:=[op(DataList),Cur]:
		fi:
		i:=i+1:
	fi:
od:

if verbose=1 or verbose=true then
	print(`Computation time:`):
	return(time()-StartTime):
fi:

DataList:
end:

###############################################
###############################################

Round:=proc(a,d:=0)
if d<=0 then return evalf(round(a)) fi:
evalf(round(a*10^d)*10^(-d)):
end:

KeepDigits:=proc(a,b) local Diff,N:
if a=b then return [a,infinity]: fi:

N:=0:
while Round(a,N)=Round(b,N) do
	N:=N+1:
od:
[Round(a,N-1),N-1]:
end:




###############################################
###############################################


FellerBook2DSafe:=proc(S,M,K1,K2,K3,K4,verbose:=1) local P,p,i,j,Div,Moms,t,k,m,X,Y,N,StartTime,
	Pow,C1,C2,C,X1,N1,X2,N2,DataList, Cur:

StartTime:=time():

if verbose=1 or verbose=true then

	Div:=`------------------------------------`:

	print(`WALKS IN THE PLANE FROM (0,0) TO (N,N)`): 
	print(``):
	print(Div):
	print(`This paper is concerned with walks from (0,0) to (N,N) where`):
	print(`steps from some subset of`,S,`are allowed.`):
	print(``):
	print(`Given a nonempty subset E of`,S,`define the generating function`):
	print(`f(N,E;t):=sum t^a(w), where w ranges over walks from (0,0)->(N,N)`):
	print(`with steps taken from E, and a is the statistic which counts`):
	print(`the number of times a walk visits {y<x} or visits {y=x} coming from {y<x}.`):
	print(``):
	print(`For each E, we generate f(N,E;t) for N=`,K1,`.. `,K2,`and `,K3,`.. `,K4):
	print(`and use this data to conjecture how moments 1 through`,M):
	print(`behave as functions of N.`):
	print(Div):
	print(``):
else
	DataList:=[]:
fi:

P:=choose(S) minus {{}}:
i:=1:
X1:=[seq(j,j=K1..K2)]:
N1:=nops(X1):

X2:=[seq(j,j=K3..K4)]:
N2:=nops(X2):

for p in P do
	Moms:=[seq(evalf(IDN(evalf( F2G([k,k],t,p) ),t,M)),k=K1..K2)]:
	for m from 1 to M do
		Y:=[seq(Moms[j][m],j=1..N1)]:
		Pow:=PowerAnsatz(m):#GuessPower(X,Y):
		C1[m]:=Aitken1(seq(Y[j]/X1[j]^Pow,j=N1-2..N1)):
	od:
	Moms:=[seq(evalf(IDN(evalf( F2G([k,k],t,p) ),t,M)),k=K3..K4)]:
	for m from 1 to M do
		Y:=[seq(Moms[j][m],j=1..N2)]:
		Pow:=PowerAnsatz(m):#GuessPower(X,Y):
		C2[m]:=Aitken1(seq(Y[j]/X2[j]^Pow,j=N2-2..N2)):
	od:
	C:=[seq(KeepDigits(C1[j],C2[j]),j=1..M)]:
	
	if {seq(C[j][1],j=1..M)}<>{0.} then
		if verbose=1 or verbose=true then
			print(`Theorem`,i,`: `):
			print(`When steps from`,p,`are allowed:`):
		else
			Cur:=[p]:
		fi:

		for m from 1 to M do
			Pow:=PowerAnsatz(m):
			if C[m][2]<>infinity then
				interface(displayprecision=C[m][2]):
			else
				interface(displayprecision=10):
			fi:
			if verbose=1 or verbose=true then
				if Pow=0 then
					print(`Moment`,m,`behaves like`,C[m][1]):
					 # `with agreement to`,C[m][2],`decimal places.`):
				elif Pow=1 then
					print(`Moment`,m,`behaves like`,C[m][1],`n`):
					  #`with agreement to`,C[m][2],`decimal places.`):
				else
					print(`Moment`,m,`behaves like`,C[m][1],`n^`,Pow):
					  #`with agreement to`,C[m][2],`decimal places.`):
				fi:
			else
				if C[m][2]=0 then Cur:=[op(Cur),infinity]:
				else Cur:=[op(Cur),C[m][1]]: fi:
				
			fi:
		
		od:
		if verbose=1 or verbose=true then
			print(Div):
		else
			DataList:=[op(DataList),Cur]:
		fi:
			
		i:=i+1:
	fi:
od:



if verbose=0 or verbose=false then
	return DataList:
fi:
print(`Computation time:`):
time()-StartTime:

end:


################################################
#Exact computations of moments
################################################

ChungFellerGF:=proc(t,k,pgf:=0,center:=0) local cfgf,cfgf1:
cfgf:=(k+1)*factorial(2*k)*((t^2)^k*t^2-1)/((t^2-1)*factorial(k+1)^2):
cfgf1:=cfgf:

if pgf=1 or pgf=true then
	cfgf1:=cfgf1/limit(cfgf,t=1):
fi:

if center=1 or center=true then
	cfgf1:=cfgf1/t^k:
fi:

simplify(cfgf1):

end:


ChungFellerMoment:=proc(m,k,standard:=1,centered:=1) local f,M,t,tD:

f:=ChungFellerGF(t,k,1,centered):

tD:=ff->t*diff(ff,t):

M:=limit((tD@@m)(f),t=1):

if standard=1 or standard=true then
	return simplify(M/ChungFellerMoment(2,k,0,1)^(m/2)):
fi:
simplify(M):
end:


ForwardKingGF:=proc(z,t):
2/(z*t+sqrt(t^2*z^2-4*t^2*z-2*t*z+1)+sqrt(z^2-6*z+1)-z):
end: